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Nomenclature  for  the  Report 

AP ,  AN ,  As,  Aw,  AE ,  AF ,  AB   Coefficients  of  the  finite  difference  equation 

m^ 
g  Acceleration  due  to  gravity,  — 

C  Scaling  factor  in  the  convection-diffusion  equation 

G  Non-dimensional  mass  flux  at  the  boundaries  of  a 

cell 

H  Height  of  the  cavity,  m 

i  x  subscript 

j  y  subscript 

k  z  subscript 

Nu  Nusselt  number 

P  Non-dimensional  pressure 

Pr  Prandd  number 

W 

q  Volumetric  heat  generation,  — 

Ra  Product  of  Rayleigh  and  Prandtl  number 

S  Source  term  in  difference  equation 

t  Non-dimensional  time 

T  Non-dimensional  temperature 

Tc  Cold  wall  temperature,  °C 

TH  Hot  wall  temperature,    °C 

AT  Temperature  difference,  Th  -  Tc,  °C 

U  Non-dimensional  x-direction  velocity 

V  Non-dimensional  y-direction  velocity 

W  Non-dimensional  z-direction  velocity 

x  Non-dimensional  horizontal  spatial  coordinate 

y  Non-dimensional  vertical  spatial  coordinate 

z  Non-dimensional  spatial  coordinate  in  the  direction 

of  depth 


Greek  Symbols 

a  Thermal  diffusivity, 


m2 
s 


(3 

K 

Ax,  Ay,Az 

v 


Coefficient  of  volume  expansion,  77- 

Generic  non-dimensional  field  variable 

W 
Thermal  Conductivity,  — ^ 

Non-dimensional  mesh  or  cell  size 

Dynamic  viscosity,  — 

kg 
Density,  -g 


Superscripts 


n 


Estimated  quantity 

Correction  to  estimated  quantity 

Present  time  level 


Subscripts 

c 

p,n,s,e,w,f,b, 

ee,ww,ss,nn, 

bb,ff 

P,  N,  S,  E,  W,F,B, 

EE,WW,SS,NN, 

BB,FF 


Refers  to  the  chip 

Designation  of  control  volume  boundaries 


Node  designation  of  basic  grid 


INTRODUCTION 


1 . 1       General  Remarks 

Over  the  years,  computational  fluid  dynamics  has  emerged  as  a  powerful  tool  to 
solve  engineering  and  scientific  problems  related  to  fluid  flow.  This  ranges  from  the 
familiar  such  as  solving  for  flow  past  a  supersonic  airplane;  flow  in  an  internal  combustion 
engine;  flow  in  an  oil  reservoir;  atmospheric  flows  and  flow  in  turbo-machines  to  the 
esoteric  such  as  geophysical  and  astrophysical  fluid  dynamics. 

The  increase  in  computing  power  which  roughly  doubles  every  two  years  has 
contributed  immensely  to  the  feasibility  of  solving  realistic  engineering  problems.  More 
stunning  advances  in  computing  resources  are  expected  in  the  near  future.  For  instance, 
massively  parallel  architecture  promises  to  create  a  prototypical  teraflop  machine  (which  is 
capable  of  1012  floating  point  operations  per  second)  by  1995  (Deng  et.  al,  1992).  This  can 
well  be  appreciated  considering  the  fact  that  three-dimensional  unsteady  problems  of 
moderate  complexity  can  be  simulated  rather  effectively  using  an  RS/6000  workstation 
which  is  capable  of  a  16  megaflops(K)6)  performance  for  a  CFD  code. 

CFD  has  gradually  evolved  to  become  a  cost-effective  alternative  to  experimentation 
with  respect  to  engineering  design  in  many  cases.  For  instance  in  the  aerospace  industry, 
the  lead  time  in  design  and  development  of  aircrafts  has  been  considerably  reduced 
(Fletcher,  1988),  thanks  to  CFD  replacing  time-consuming  and  costly  testing  procedures  in 
many  of  the  steps.  This  cost-effectiveness  is  expected  to  improve  relentlessly,  as 
computing  costs  decrease. 

Nevertheless,  CFD  must  be  used  in  the  design  process  with  utmost  care.  In  the  first 
place,  CFD  can  never  claim  to  replace  experiments;  at  least  not  yet.  Experiments  form  a 
critical  component  of  the  overall  process.  It  is  imperative  that  the  CFD  code  be  validated  by 
experiments  for  the  particular  case.  There  are  several  reasons  for  that.  The  validity  of  the 
boundary  conditions,  the  refinement  of  the  grid  or  the  validity  of  certain  simplifying 
assumptions  can  only  be  checked  and  confirmed  by  experiments.  There  might  even  be  bugs 
which  need  to  be  fixed.  For  computer  codes  which  include  turbulence  modelling  it  is 
sometimes  required  to  check  whether  the  emperically  determined  constants  are  valid  for  the 
given  case. 


The  thermal  management  of  electronic  components  (otherwise  known  as  electronic 
cooling)  can  also  benefit  from  these  advances  in  computing  power.  The  design  process  has 
in  the  past  relied  heavily  on  empirical  testing.  Convective  heat  transfer  modelling  was 
expensive  due  to  the  complex  mathematical  equations  involved.  The  heat  transfer 
coefficients  were  determined  from  experiments  and  were  then  used  in  simpler  conduction 
calculations  to  determine  the  maximum  chip  temperature  and  other  quantities  of  interest  to 
an  engineer. 

With  advances  in  Computational  Fluid  Dynamics  and  computing  power,  the 
convection  and  conduction  problems  can  be  solved  simultaneously  in  a  coupled  manner. 
Computational  fluid  flow  and  heat  transfer  is  increasingly  being  used  as  part  of  this  design 
process.  A  number  of  general-purpose  commercial  code  available  (such  as  PHOENICS, 
FLUENT,  FIDAP  etc)  can  be  used  to  solve  applications  in  electronic  cooling.  There  are 
several  drawbacks  though.  The  codes  are  too  general  (and  very  expensive!).  A 
considerable  amount  of  effort  must  be  expended  in  order  to  tailor  them  to  the  specific 
application.  More  important,  the  source  code  is  almost  never  available.  The  option  of 
modifying  the  code  does  not  arise. 

The  present  code  AMPHIB  attempts  to  redress  these  disadvantages.  It  is  a  general- 
purpose  code  that  is  exclusively  designed  to  solve  problems  in  electronic  cooling. 
AMPHIB  is  also  fully  validated  by  experiments.  It  can  therefore  be  used  directly  in  such 
applications  with  minimal  modifications  if  at  all.  To  the  best  of  our  knowledge,  this  is  the 
first  software  package  that  specifically  solves  these  class  of  problems.  Since  the  source 
code  is  readily  available,  the  code  can  be  modified  to  solve  problems  of  much  greater 
complexity  such  as  those  involving  turbulence  modelling,  mass  transfer  etc. 

1 . 2       Purpose  of  the  Code 

The  primary  purpose  of  the  code  is  as  a  design  tool  for  applications  in  electronic 
cooling.  As  noted  earlier,  experiments  form  an  integral  part  of  this  procedure.  The 
usefulness  of  the  code  is  best  illustrated  with  examples. 

Consider  a  hypothetical  configuration  shown  in  figure  1.1.  A  power  dissipating 
chip  is  protruding  from  a  vertical  substrate.  It  is  being  cooled  by  a  dielectric  fluid  in  direct 
contact  with  the  chip  surface.  The  dielectric  fluid  is  hermetically  sealed  in  an  enclosure 
which  contains  the  chip  and  substrate.  Heat  exchangers  at  the  top  and  the  bottom  wall 
extract  heat  from  the  liquid.  This  configuration  is  what  is  commonly  known  as  cooling  by 
liquid  immersion. 


To  check  design,  experiments  may  have  been  carried  out  to  monitor  the  chip 
temperature  for  a  case  that  is  considered  optimum.  The  chip  symmetrically  placed  at  the 
center  of  the  substrate  would  seem  to  be  a  reasonably  optimal  design.  Chip  temperatures 
are  recorded  for  different  power  levels  in  the  chip.  However,  some  nagging  doubts  still 
remain.  Is  this  really  the  optimal  configuration?  What  if  the  chip  was  placed  in  an  eccentric 
position?  Would  the  chip  temperature  be  reduced?  Also,  would  it  matter  if  the  substrate  was 
a  different  material  of  higher  conductivity  e.g,  ceramic  instead  of  plexiglas.  What  would  be 
the  effect  of  increasing  the  spacing  between  the  substrate  and  the  vertical  wall  opposite  to 

it? 

These  are  legitimate  questions.  Unfortunately  it  is  very  time  consuming  to  respond 
to  these  questions  experimentally.  This  is  where  CFD  codes  such  as  AMPHIB  fit  in.  As  a 
first  step,  the  code  is  used  to  validate  the  experiments.  The  chip  temperatures  predicted  by 
the  program  must  be  made  to  match  those  that  were  experimentally  determined  within  an 
acceptable  bound. 

Appropriate  models  of  the  boundary  conditions  must  be  introduced  in  the  program. 
To  have  a  well  defined  problem  the  geometrical  and  physical  parameters  of  the  setup  must 
be  included  in  the  code.  This  includes  the  complete  dimensions  of  the  enclosure,  the 
dimensions  of  the  chip,  the  thermo-physical  properties  of  the  chip,  substrate  and  fluid.  The 
power  dissipated  by  the  chip  and  the  temperature  of  the  heat  sinks  (  the  top  and  bottom 
wall)  and  the  location  of  the  chip  must  also  be  known. 

The  present  code  (AMPHIB)  can  handle  the  problem  so  described  without 
modifications.  If  the  code  predicts  reasonably  well,  it  is  quite  straightforward  to  change 
the  parameters  to  predict  chip  temperatures  for  a  different  geometry  or  different  materials. 
The  exact  details  of  the  code  and  implementation  are  given  in  a  different  chapter. 

Consider  another  configuration  of  a  greater  complexity  shown  in  figure  1.2.  A 
three  by  three  array  of  chip  protrusions  is  shown.  Experiments  were  performed  with  a 
certain  spacing.  All  the  chip  dimensions  are  the  same.  The  chip  is  so  configured  that  the 
largest  dimension  is  aligned  to  the  vertical.  Experiments  were  performed  for  a  certain 
spacing  between  the  substrate  and  the  opposite  wall.  Unfortunately,  due  to  constraints  in 
the  instrumentation,  experimentation  was  possible  only  for  spacings  that  were  larger  than 
the  design  constraint.  Spacing  between  chips  and  wall  are  typically  very  small  since  the 
chip  density  in  real  applications  are  quite  high. 

The  numerical  code  can  be  very  usefully  applied  to  such  a  situation.  Firstly,  the 
code  validates  the  experimental  results  for  the  large  spacing.  Later,  numerical  computations 
are  done  for  the  smaller  design  spacing.  It  is  also  possible  to  check  computationally 
whether  a  slightly  different  alignment  of  chips  is  a  good  idea.  For  instance,  if  the  largest 
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dimension  of  the  identical  chips  was  aligned  to  the  horizontal  rather  than  the  vertical 
direction.  This  is  a  rather  trivial  thing  to  do  computationally  but  a  rather  cumbersome  thing 
to  do  experimentally.  In  the  latter  case,  a  whole  new  experimental  setup  has  to  be 
constructed.  The  optimal  approach  is  to  judiciously  combine  numerics  with  experimentation 
in  the  design  process.  The  configuration  shown  in  figure  1.2  can  be  handled  by  AMPHIB 
with  no  changes.  In  fact,  the  program  can  handle  an  n  x  m  array  of  chips  on  a  substrate 
where  n  and  m  are  whole  numbers.  Theoretically,  there  are  no  limits  on  the  number  of 
chips  in  the  computations.  Nevertheless,  a  5  by  5  chip  array  would  itself  require  107  bytes 
of  memory  for  reasonably  accurate  results  in  three-dimensional  computations.  Thus,  there 
are  practical  limits  on  the  total  number  that  can  be  handled  by  the  code. 

There  are  other  related  configurations  which  can  be  tackled  by  the  program  with 
minor  modifications.  One  case,  which  is  important  from  the  engineering  point  of  view  is 
shown  in  figure  1.3.  The  array  of  chips  shown  is  cooled  by  a  fluid  that  is  circulated  by  a 
pump.  This  is  a  mixed  convection  problem.  An  important  question  is  should  the  fluid 
circulate  from  the  top  to  the  bottom  or  the  other  way  round.  The  answer  may  not  be  simple. 
It  could  depend  on  the  power  dissipation  rate,  the  properties  as  well  as  the  geometry  of  the 
configuration.  The  code  can  be  used.  The  modifications  in  the  code  that  are  necessary  to 
solve  such  a  problem  are  discussed  in  a  later  chapter. 

A  more  general  configuration  is  shown  in  figure  1.4.  The  whole  chip  assembly  and 
the  enclosure  is  at  angle  with  the  vertical  and  fluid  is  circulated  through  it.  This  problem  can 
be  solved  with  surprisingly  little  changes  with  AMPHIB. 

1 . 3       Overview  of  the  Report 

In  chapter  two,  the  algorithm  of  the  code  is  presented.  The  description  is  brief. 
Interested  readers  can  look  into  the  references  provided  in  that  chapter  for  additional  details. 
A  good  readable  account  of  the  control  volume  method  applied  to  heat  transfer  and  fluid 
flow  problems  is  given  in  Patankar  (1980).  Readers  interested  in  the  code  purely  as  a  user 
are  advised  to  go  to  chapter  four  directly  which  describes  the  inputs  and  subroutines  of 
AMPHIB. 

Chapter  three  describes  three  example  problems  that  were  looked  into.  One  of  the 
example  deals  with  convection  and  conduction  in  an  enclosure  with  a  three  by  three  array  of 
chips  embedded  in  a  substrate.  Another  example  is  convection  in  an  enclosure  with  no 
substrate  or  chip.  Both  these  examples  are  compared  with  experiments.  That  is  precisely 
why  they  were  chosen.  These  two  examples  underscore  the  essential  correctness  and 
accuracy  of  the  code.  The  third  example  is  not  related  to  any  experiments.  It  is  a 


benchmark.  It  is  useful  to  have  one,  especially  when  changes  are  made  in  the  code.  One 
can  always  fall  back  on  this  case  for  debugging  and  checking  purposes. 

Chapter  four  provides  a  detailed  and  complete  account  of  the  inputs  to  the  program 
(all  of  them  initialized  in  BLOCK  DATA)  and  the  subroutines.  Appendix  A  details  the 
implementation  of  the  boundary  conditions.  This  chapter  is  important  for  someone  who 
intends  to  modify  the  code  and  thus  extend  the  scope  of  AMPHIB.  The  present  version  of 
the  code  has  a  fixed  type  of  boundary  condition.  There  are  no  options  or  flags  in  the  code 
which  would  enable  one  to  chose  an  arbitrary  set  of  boundary  conditions  for  each  wall. 
Appendix  A  has  easy  to  follow  instructions  on  how  to  change  the  code  in  order  to 
accommodate  a  more  general  set  of  boundary  conditions. 

Appendix  B  provides  explicit  algebraic  expressions  for  the  difference  equations  and 
numerical  boundary  conditions  that  are  implemented  in  the  code.  Appendix  C  gives  a  listing 
for  AMPHIB. 
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Heat  Exchanger 


Heat  exchanger 


Fluid-filled 
enclosure 


Figure  1.1 

In  this  configuration,  a  single  chip  embedded  in  a  substrate  is 
being  cooled  by  a  dielectric  fluid.  The  fluid  is  hermetically  sealed  in 
a  chamber.  The  heated  fluid  is  cooled  by  heat  exchangers  at  the  top 
and  bottom  walls.  The  first  step  is  to  mathematically  model  the 
physical  problem.  For  example,  the  top  and  bottom  wall  can  be 
modelled  as  isothermal  heat  sinks.  The  computer  code  can  then  be  used 
to  calculate  the  chip  temperatures  for  different  geometrical  parameters 
and  transport  properties. 


1  1 


Heat  Exchanger 


Substrate 


Chip  Locations 


Dielectric  Fluid 
filled  enclosure 


Heat  Exchanger 


Figure  1.2 

This  is  a  chip- substrate  configuration  of  a  slightly  greater 
complexity.  There  are  nine  chips  in  a  three  by  three  array 
that  is  being  cooled  by  a  dielectric  fluid  in  an  enclosure. 
This  problem  can  be  solved  by  AMPHIB  without  changes 
provided  that  the  chip  locations  are  along  a  rectangular  grid. 
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Fluid  inlet 


Fluid  Outlet 


Figure  1.3 


This  is  a  configuration  similar  to  figure  1.2,  except 

that  cooling  is  now  aided  by  fluid  that  is  being  circulated 

by  a  pump.  This  problem  can  be  solved  by  AMPHIB 

only  after  modifications  are  being  made  in  the  code  which  is 

described  at  length  in  a  later  chapter. 
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Fluid  inlet 


Substrate 


Chip 


ssssssssssssssgssssgss 


} 


t 


Fluid  Outlet 


Figure  1.4 


The  direction 
of  the  gravity 
vector 


This  is  a  chip- substrate  configuration  that  is  being 
cooled  in  a  mixed-convection  mode.  However,  as  in 
many  real-life  situations,  the  whole  assembly  is  at 
an  angle  with  the  vertical.  Nevertheless,  this  problem 
can  be  solved  by  AMPHD3  with  minimal  changes  in  the 
program. 
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OUTLINE  OF  ALGORITHM 


2. 1       The  Scope  of  the  Program 

The  program  is  a  general  purpose  finite  difference  FORTRAN-77  code  designed  to 
solve  unsteady  and  steady  three-dimensional  heat  transfer  and  fluid  flow  problems  in 
rectangular  cartesian  coordinates.  The  code  is  designed  to  effectively  handle  conjugate  heat 
transfer  situations,  i.e.,  when  conduction  and  convection  is  solved  simultaneously. 
However,  pure  convection  or  conduction  problems  can  also  be  solved. 

The  program  uses  the  control  volume  method  (Patankar,  1980)  to  discretize  the 
governing  Navier-Stokes,  continuity  and  the  energy  equations  in  the  primitive  variable 
formulation.  The  SIMPLEX  algorithm  (Van  Doormal  and  Raithby,  1985)  is  used  to 
determine  the  five  field  variables  (the  three  velocity  components,  temperature  and 
pressure).  It  is  essentially  a  more  implicit  version  of  SIMPLE  (Patankar,  1980).  Due  to  the 
more  implicit  nature  of  the  algorithm,  larger  time  steps  can  be  used  without  encountering 
numerical  instability.  This  leads  to  considerable  saving  of  computer  resources  for  steady 
state  computations  since  a  smaller  number  of  time  steps  would  be  required.  The 
computational  work  per  time  step  would  be  increased  somewhat,  but  is  more  than 
compensated  by  the  decrease  in  the  total  number  of  time  steps  required  to  reach  a  converged 
solution. 

The  algorithm  is  for  three-dimensional  and  unsteady  heat  transfer  and  fluid  flow 
problems  that  marches  in  time.  Time  stepping  is  performed  by  a  first  order  Back  ward- Euler 
scheme.  The  algorithm  is  therefore  implicit  in  time.  The  grids  are  staggered,  and  the 
QUICK  scheme  is  used  for  interpolating  the  convective  terms  (Leonard,  1983).  Staggered 
grid  means  that  grid  points  for  the  velocities  (the  so  called  vector  grid)  are  different  from 
those  of  pressure  and  temperature  grid  points  (the  scalar  grid).  More  will  be  said  of  it  later. 
Generally  speaking,  the  staggered  grid  leads  to  more  accurate  results  since  the  need  to 
interpolate  reduces  somewhat.  It  is  now  an  accepted  norm  for  all  CFD  calculations  with 
finite  differences  with  the  primitive  variable  formulation. 

QUICK  (quadratic  interpolation  for  convective  kinetics)  is  a  third  order  polynomial 
interpolation  scheme  for  the  convective  terms.  It  is  almost  as  stable  as  the  first  order 
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upwind  scheme  and  vastly  more  accurate.  It  does  have  the  twin  advantages  of  being  robust 
as  well  as  accurate  compared  to  the  more  conventional  techniques  such  as  upwinding  or 
central  differencing  which  are  the  mainstay  of  most  if  not  all  commercially  available  codes 
to  solve  convection  problems.  It  is  a  well  known  fact  (Patankar,  1988)  that  although 
QUICK  is  very  accurate  it  is  not  as  robust  as  the  upwind  scheme.  There  are  times  when  the 
QUICK  scheme  might  not  converge  due  to  insufficient  number  of  grid  points  in  regions  of 
high  velocities.  The  code  therefore  has  the  option  of  using  the  upwind  scheme  when  such  a 
situation  arises.  In  particular,  the  code  uses  a  scheme  which  is  a  weighted  average  of 
QUICK  and  upwind.  Depending  on  the  case  one  can  have  a  range  of  schemes  ranging  from 
pure  upwind  to  pure  QUICK. 

Except  for  the  viscosity,  the  fluid  is  assumed  to  have  constant  transport  properties 
and  is  incompressible  except  for  density  changes  to  account  for  buoyancy.  Viscous 
dissipation  and  pressure  work  are  neglected.  These  assumptions  are  justified  for  most 
applications.  The  harmonic  mean  formulation  for  the  interfacial  diffusivities  is  used  to 
effectively  handle  sharp  discontinuities  in  property  values  in  the  computational  domain  that 
arise  in  conjugate  problems  (Patankar,  1980).  Such  is  the  case  in  electronic  cooling 
problems  since  the  thermo-physical  properties  vary  drastically  between  the  fluid  and  the 
solid  regions.  A  substrate  and  package  are  examples  of  solid  regions  where  properties 
might  differ  sharply. 


2.2       Algorithmic  Details 

The  code  uses  the  Control  Volume  Method  to  solve  five  coupled  non-linear  partial 
differential  equations  (governing  equations).  All  the  equations  can  be  cast  in  the  form  of  a 
three-dimensional  convection-diffusion  equation  of  the  following  type: 

The  different  governing  equations  can  be  recovered  by  setting  different  values  for  the 
canonical  variable  0  and  the  constant  C  as  well  as  the  source  term  S  which  are  tabulated 
below. 


Equation  (() 
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ap 

x-momentum  U  Pr  "  ^x 

dP 
y-momentum  V  Pr  -  ^     +  T  Pr  Ra 

3P 
z-momentum  W  Pr  -  ^z 

Energy  T  1  0 

Continuity  1  0  0 

The  governing  equations  are  a  result  of  a  very  specific  set  of  non-dimensionalization  of  the 
dependent  and  independent  variables.  The  x,  y  and  z  coordinates  are  non-dimensionalized 

a 

by  H,  the  enclosure  height.  The  velocities  are  scaled  by  pj ,  the  time  is  non-dimensionalized 

by  — ,  the  pressure  difference  by  ^p-;  and  the  temperature  difference  by  ^.  q  is  the 

volumetric  rate  of  heat  generation  per  chip,  a  is  the  thermal  diffusivity  of  the  fluid,  p  is  the 
fluid  density  and  the  the  conductivity  of  the  chip  is  kc.  For  a  different  set  of  scales  the 
governing  equations  will  be  of  the  same  form  as  equation  (1)  but  with  different  values  of  C 

and  S. 

A  grid  is  imposed  in  the  physical  domain  and  the  governing  equations  are 
discretized  on  the  grid.  For  the  control  volume  method,  the  governing  equations  are  first 
integrated  around  a  volume  about  the  grid  point  known  as  the  control  volume.  The  problem 
is  now  reduced  to  that  of  calculating  the  values  of  the  dependent  variables  at  the  grid 
locations.  The  generic  convection-diffusion  equation  (equation  1)  is  reduced  to  an  algebraic 
equation  of  the  following  form: 

Ap<t>p  =  A^  +  Aw^  +  An<j>n  +  As(|>s  +  AE<|>E  +  Ab<j>b  +  Ap<j)F  +   S        (2.2) 

The  scalar  variables  (pressure  and  temperature)  are  in  the  non-staggered  grid.  The  grid 
consists  of  the  grid  points  and  control  volume  surrounding  it.  In  figure  2. 1  the  shaded 
region  indicates  the  non-staggered  grid.  The  grid  points  are  indicated  by  the  circular  dots. 
The  subscripts  refer  to  the  location  with  respect  to  the  grid  point  as  shown  in  the  same 
figure. 

1.  The  grid  point  in  question  is  referred  to  by  the  subscript  P. 

2.  The  grid  point  to  the  right  of  P  is  subscripted  by  E  (east).  Similarly  the  grid  point  to  the 
right  of  E  is  the  EE  grid  point  as  indicated  in  figure  2.2. 

3.  The  grid  point  to  the  left  of  P  is  the  W  (west)  point. 
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4.  The  grid  point  above  P  is  N  (north). 

5.  The  grid  point  below  P  is  S  (south). 

6.  The  B  and  F  grid  points  are  not  shown  in  the  figure.  The  B  (back)  grid  point  is  located 
above  the  paper  immediately  above  P.  This  is  consistent  with  a  right-handed  cartesian 
system.  Similarly,  the  F  (front)  grid  point  is  located  under  the  plane  of  the  paper. 

7.  The  grid  points  WW,  NN,  SS,  BB  and  FF  are  similarly  defined. 

The  grids  are  chosen  such  that  the  walls  of  the  enclosure  coincide  with  the  control 
volume  faces  of  the  non-staggered  grid.  The  staggered  grids  are  for  the  vector  variables 
(the  three  velocities  u,  v  and  w).  The  staggered  grids  are  shifted  by  half  a  control  volume 
towards  their  respective  negative  directions.  In  figure  2.1  the  grid  point  for  the  u- velocity 
(velocity  in  the  x-direction)  are  shown  as  triangular  dots.  The  control  volume  surrounding 
the  u-velocity  grid  is  highlighted  by  bold  dashed  lines. 

The  u-velocity  grid  points  are  thus  located  at  the  east  and  west  faces  of  the  non- 
staggered  control  volume.  The  staggered  grid  points  for  the  v-velocity  are  portrayed  as 
square  dots.  The  staggered  v-velocity  control  volume  is  enclosed  by  bold  dashed  lines  just 
as  the  other  one.  The  v-velocity  grid  points  coincide  with  the  north  and  south  face  of  the 
scalar  grid  as  shown  in  the  figure. 

The  coupled  fluid  flow  and  heat  transfer  problem  is  solved  using  the  SIMPLEX 
algorithm.  Details  are  given  in  Van  Doormal  and  Raithby  (1985).  Briefly,  equations  are 
solved  for  the  temperature,  the  three  velocities,  the  pressure  correction  and  the  coupling 
constants  for  the  three  momentum  equations.  The  flow  chart  is  shown  in  figure  2.2.  The 
temperature,  the  velocities,  the  coupling  constants,  and  the  pressure  corrections.  The  most 
updated  field  variables  are  used  at  each  stage,  e.g.,  for  solving  the  v-velocity  the  updated  u- 
velocity  and  temperature  are  used  and  so  on. 
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Figure  2.1 

The  staggered  grid,which  is  one  of  the  cornerstones  of  the  control  volume  method 

is  shown.  The  non-staggered  control  volume  is  bounded  by  plain  lines.  The  dotted 

lines  pass  through  the  grid  points.  The  non-staggered  grid  point  which  is 

the  location  of  the  temperature  and  pressure  is  represented  by  circular  dots.  The 

u-velocity  grid  points  are  represented  by  triangular  dots.  Similarly,  the  v-velocity 

grid  points  are  represented  by  square  dots.  The  grid  in  the  third  direction  is 

analogous. 


19 


START 


I 


t  =  0 


Specify  initial  field  of  the 
temperature,  velocities 
and  pressure 


1 


t  =  t  +  dt 


I 


No 


Assemble  the  energy 
equation  and  solve  for 
temperature  (T) 


I 


No 


Assemble  x-mom  equation 
Solve  U  (x-velocity) 
Solve  DU  (x-mom  coupling 

rnrKtant^ 


1SL 


Assemble  pressure  correction 
eqn  and  solve.  Correct 
pressure,  velocities  and 
calculate  residual  mass. 


Assemble  y-mom  equation 
Solve  V  (y-velocity) 
Solve  DV  (y-mom  coupling 
constant) 
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Assemble  z-mom  equation 
Solve  W  (z- velocity) 
Solve  DW  (z-mom  coupling 
constant) 


Figure  2.2 

The  coupled  heat  transfer  and  fluid  flow  problem  is  solved  by  an  iterative  algorithm 
known  as  SIMPLEX.  Iterations  are  necessary,  since  the  problem  is  non-linear.  The 
algebraic  equations  are  linearized  at  every  iteration  step,  and  are  approximate.  Each 
iteration  provides  an  improvement  on  this  approximation.  After  many  iterations,  a 
converged  and  accurate  solution  of  the  non-linear  problem  is  obtained.  In  contrast, 
if  only  conduction  is  solved  for,  no  iterations  are  needed  since  the  equations  of 
conduction  are  linear. 


20 


EXAMPLES 


3.1  General  Remarks 

Some  specific  examples  are  discussed  in  this  chapter.  This  serves  as  a  validation  for 
the  code  for  the  given  applications,  namely  cooling  of  a  three  by  three  array  of  chips 
mounted  on  a  substrate  as  well  as  convection  without  chip  or  substrate. 

3.2  Convection  in  the  absence  of  Substrate  and  Chip 

The  boundary  conditions  are  the  following: 

1.  Top  wall  cooled  and  bottom  wall  heated.  THOT  =  0.5,  TCOOL  =  -  0.5. 

2.  All  side  walls  insulated. 

The  parameters  as  initialized  in  BLOCK  DATA  are  the  following: 

DATA  NIP2,  NIP1,  NI,NIM  1/37,36,35,34/ 

DATA  NJP2,N JP1  ,NJ,NJM  1/33,32,3 1 ,30/ 

DATA  NKP2,NKPl,NK,NKMl/27,26,25,24/ 

DATA  NNMAX,NMAX,IMAX,ITMAX,KRUN,NPRINT/29952,5000, 10,5, 1,100/ 

DATA  SMALL,EPS,SORMAX/l .0E20, 1  .OE-8,2.0/ 

DATA  ISUB,ICHIP,JCHIP,KCHIP,NCHP,ICHOICE/2,2,2,2,0,0/ 

DATA  XBR,YBR,ZBR/0.40,0. 13,0.30/ 

DATA  H,WTH,BTH/12.3,14.6,25.0/ 

DATA  HCHIP,WCHIP,BCHIP,BSUB/24.0,8.0,6.0,19.5/ 

DATA  ALC,ALS,THOT,TCOOL,TAVG/2823.0,3.5,0.5,-  0.5,19./ 

DATA  RHS,RHC/1.20,0.262/ 

DATA  QQQ,DTIME,RA  ,PR,XPER,ROLL/1.5,1.0E-5,1.0E7,130.,2.0,2.0/ 

DATA  IUNFRM/0/ 

DATA  QUICK/1.0/ 

DATA  NJCHIP,NKCHIP/3,3/ 

DATA  YCHIP,ZCHIP/38.0,76.0,114.0,50.8,101.6,152.4/ 
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Note  that  the  number  of  prescribed  parameters  are  in  excess  of  what  is  required.  For 
example,  when  NCHP  is  set  to  0  the  chip  parameters  such  as  ISUB,  HCHIP  etc  are  not 
required  and  can  be  set  quite  arbitrarily. 

In  figure  3.2,  the  numerical  results  are  compared  with  experiments.  Both  the 
numerical  and  experimental  results  are  for  a  natural  convection  problem  with  no  substrate 
or  chip  in  a  three-dimensional  rectangular  enclosure.  The  heated  bottom  wall  provides  the 
driving  force  for  the  flow.  Figure  3.1  provides  a  schematic  diagram  of  the  enclosure 
geometry.  Figure  3.2  shows  the  velocity  field  across  a  vertical  mid- section  as  shown  in 
figure  3.1.  The  results  from  four  different  cases  are  compared  with  the  experiments  due  to 
Arroyo  and  Saviron  (1992). 

The  velocity  vectors  on  the  right  side  are  the  experiments  and  on  the  left  are  the 
numerical  computations.  Each  pair  corresponds  to  a  particular  Rayleigh  number.  The 
BLOCK  DATA  has  been  initialized  for  a  Rayleigh  number  of  44,744.  The  Rayleigh 
numbers  for  the  five  cases  starting  from  the  top  are  the  following:  6319.0,  11101.0, 
22884.0,  44744.0  and  66604.0.  For  high  Rayleigh  number  computations,  it  is  generally 
sound  practice  to  begin  at  a  lower  Rayleigh  number  and  work  ones  way  upwards.  In  other 
words,  use  a  lower  Rayleigh  number  as  an  initial  condition  to  start  a  higher  Rayleigh 
number  calculation.  Also,  note  that  frequent  restarts  are  necessary  even  for  a  particular 
Rayleigh  number  since  steady-state  may  not  be  achieved  in  a  single  run. 

Figure  3.3  compares  the  computations  (on  the  left)  with  the  same  set  of  experiments 
for  the  first  three  cases  ( i.e.,  6319.0,  11101.0  and  22884.0).  What  is  now  compared  are 
the  same  set  of  velocities.  However,  the  data  is  displayed  in  terms  of  pathlines.  Since, 
pathlines  can  be  observed  by  very  standard  flow  visualization  techniques,  it  is  a  very  useful 
thing  to  compare.  The  degree  of  match  between  the  two  is  striking.  Also  note  that  the  code 
does  not  need  to  be  modified.  It  must  be  emphasized,  that  such  a  good  comparison  between 
experiments  and  numerics  are  very  rare.  It  is  not  often  that  such  a  comparison  is  even 
attempted  for  three-dimensional  flows.  Thus,  the  results  do  demonstrate  the  essential 
correctness  of  the  code. 

3.2       Convection  with  Substrate  and  Chip 

The  boundary  conditions  are  the  following: 

1.  The  top  and  bottom  wall  are  isothermal  heat  sinks.  This  case  is  the  mathematical 
modelling  of  the  case  shown  in  figure  1.2.  THOT  =  0.0,  TCOOL  =  0.0. 

2.  All  side  walls  insulated. 

The  parameters  as  initialized  in  BLOCK  DATA  are  the  following: 
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DATA  NIP2,  NIP1,  NI,NIM1/21, 20,19, 18/ 

DATA  NJP2,NJPl,NJ,NJMl/39,38,37,36/ 

DATA  NKP2,NKPl,NK,NKMl/66,65,64,63/ 

DATA  NNMAX,NMAX,IMAX,ITMAX,KRUN,NPRINT/49400,4000, 10,5, 1 ,  100/ 

DATASMALL,EPS,SORMAX/1.0E20,1.0E-8,2.0/ 

DATA  ISUB,ICfflP,JCHIP,KCHIP,NCHP,ICHOICE/4,4,4,3,l,l/ 

DATA  XBR,YBR,ZBR/0.75, 1.00, 1.00/ 

DATA  H,WTH,BTH/152.,203.2,49.5/ 

DATA  HCHIP,WCHIP,BCHIP,BSUB/24.0,8.0,6.0,19.5/ 

DATAALC,ALS,THOT,TCOOL,TAVG/3338.0,2.0,0.0,0.0,19.0/ 

DATAQQQ,DTIME,PR,RA,XPER,ROLL/1.5,2.0E-8,200.0,1.0E9,0.0,0.0/ 

DATA  IUNFRM/0/ 

DATA  QUICK/1.0/ 

DATA  NJCHIP,NKCHIP/3,3/ 

DATA  YCHIP,ZCHIP/38.0,76.0,114.0,50.8,101.6,152.4/ 

This  corresponds  to  a  case  with  a  three-by-three  array  of  chips  that  are  configured 
uniformly  and  symmetrically  on  one  of  the  vertical  sidewalls  that  includes  the  substrate  as 
well.  The  run  is  for  a  fairly  fine  grid  (20x38x65).  The  fluid  is  the  Fluorinert  liquid  FC71. 
The  viscosity  is  assumed  to  be  a  strong  function  of  temperature.  The  power  dissipation  per 
chip  is  1.5W.  The  flow  field  for  this  set  of  computations  is  highly  three-dimensional  and 
unsteady. 

The  dimensions  of  chips  and  enclosure  are  shown  in  figure  3.4.  Figure  3.5  shows 
some  of  the  preliminary  data  that  have  been  obtained  from  such  a  study.  For  a  FC75  fluid  it 
can  be  seen  in  figure  3.5(a)  that  the  heat  lost  to  the  substrate  by  conduction  from  the  chip 
decreases  with  an  increase  in  the  Rayleigh  number.  This  is  not  surprising,  since  high 
Rayleigh  number  flows  will  be  more  convection  dominated.  For  higher  Rayleigh  numbers 
an  increasing  proportion  of  the  heat  is  transferred  directly  to  the  fluid.  The  substrate  is 
made  of  plexiglas.  Things  would  of  course  be  different  if  the  substrate  was  a  good 
conductor  such  as  ceramic. 

Figure  3.5(b)  shows  the  chip  temperatures  for  the  same  fluid  and  same  geometry  as 
a  function  of  the  Rayleigh  number.  It  is  interesting  to  note  that  a  straight  line  is  obtained  in 
the  log-log  scale.  It  is  therefore  quite  likely  that  a  simple  correlation  can  be  deduced  which 
can  later  be  applied  for  design  purposes. 
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Figure  3.6  compares  the  experimentally  obtained  pathlines  due  to  Joshi  et  al.  (1990) 
with  the  computations  using  AMPHIB  for  the  case  initialized  in  BLOCK  DATA.  The 
power  level  is  1.5  Watts  per  chip  and  the  fluid  is  FC71.  The  viscosity  varies  by  more  than 
a  factor  of  10.  This  is  because  the  viscosity  is  a  strong  function  of  the  temperature  for  the 
given  fluid.  One  observes  a  relatively  stagnant  flow  near  the  bottom  and  vertical  plumes 
along  the  chips.  A  coarser  grid  would  not  suffice  for  this  problem,  since  the  thin  boundary 
layer  along  the  chips  and  the  walls  need  to  be  resolved.  The  results,  do  indicate  that 
AMPHIB  could  be  usefully  employed  as  a  design  tool.  However,  powerful  computational 
resources  are  required.  A  number  cruncher  of  the  class  of  a  RS/6000  workstation  is  a  must 
for  realistic  three-dimensional  calculations  such  as  these. 

3.4       Convection  without  Chip  or  Substrate:  Benchmark  Study 

This  example  is  set  up  as  a  benchmark.  The  boundary  conditions  are  the  following: 

1.  Top  wall  cooled  and  bottom  wall  heated.  THOT  =  0.5,  TCOOL  =  -  0.5. 

2.  All  side  walls  insulated. 

The  parameters  as  initialized  in  BLOCK  DATA  are  the  following: 

DATA  NIP2,  NIP1,  NLNIMl/13,12,11,10/ 

DATA  NJP2,NJP1,NJ,NJM1/13,12,11,10/ 

DATA  NKP2,NKP1,NK,NKM1/13,12,11,10/ 

DATANNMAX,NMAX,IMAX,ITMAX,KRUN,NPRINT/1728,500,10,5,0,100/ 

DATA  SMALL,EPS,SORMAX/l .0E20, 1 .0E-8,2.0/ 

DATA  ISUB,ICHIP,JCHIP,KCHIP,NCHP,ICHOICE/2,2,2,2,0,0/ 

DATA  XBR,YBR,ZBR/0.75,1.00,1.00/ 

DATA  H,WTH,BTH/100.,100.,100./ 

DATA  HCHIP,WCHIP,BCHIP,BSUB/24.0,8.0,6.0,19.5/ 

DATAALC,ALS,THOT,TCOOL,TAVG/3338.0,2.0,0.5,-0.5,25.0/ 

DATA  QQQ,DTIME,PR,RA,XPER,ROLL/1.5,2.0E-3,2.5,2.5E4,0.0,0.0/ 

DATA  IUNFRM/1/ 

DATA  QUICK/1.0/ 

DATA  NJCHIP,NKCHIP/3,3/ 

DATA  YCHIP,ZCHIP/38.0,76.0,1 14.0,50.8,101.6,152.4/ 

After  500  time  steps,  at  the  end  of  the  computations,  the  heat  and  mass  balance  statistics 
looked  like  the  following: 
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NT=      500  TIME=  1.0000 

ITER=  1  SOURCE=  0.001057  SORSUM  =  0.000000 

NUC  =  2.117941  NUH=2. 117934 

QCHIP  =  0.00000  QALL=  -  0.000007 

QUNS  =  -  0.001424  QCRT  =  0.000002 

QUNS  is  relative  degree  of  unsteadiness  in  the  temperature.  For  steady-state  to  be  reached, 
QUNS  should  be  about  106  of  NUC. 

This  doesn't  correspond  to  any  experimentally  observed  case.  It  is  nevertheless  a  useful 
benchmark.  This  case  doesn't  take  much  computer  time.  It  runs  for  about  221  seconds  in 
the  CONVEX  C240,  90.2  seconds  in  the  AMDAHL  5990/500,  and  about  66.2  seconds  in 
the  CRAY  X-MP/216.  If  large  scale  modifications  have  been  made  in  the  style  and  content 
of  the  code  it  will  be  useful  to  run  the  benchmark  and  check  to  see  if  the  numbers  are 
unchanged  in  the  output  and  as  such  ensures  that  the  changes  have  been  free  of  errors.  It 
must  be  emphasized  that  the  numbers  after  modifications  must  be  unchanged.  It  is  our 
experience  that  relatively  serious  bugs  cause  minor  changes  in  the  output. 

Note  that  the  inputs  ROLL  and  XPER  will  not  be  described  in  the  next  chapter. 
They  are  basically  not  relevant  to  computations  in  liquid  cooling.  XPER  must  be  set  to  0. 
and  ROLL  can  be  set  arbitrarily. 
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Top  wall  cooled  isothermally 


Gravity  Vector 


Velocity  vector  for 
subsequent  plots  are 
in  this  section. 


All  side  walls  insulated  including 
this  one. 


Bottom  wall  isothermally  heated 


Figure  3.1 


This  is  an  illustration  of  a  problem  with  a  simple  geometry  but  complicated  physics,  i.e, 
Rayleigh-Benard  Convection.  The  driving  force  is  basically  the  hot  bottom  wall.  The 
top  wall  is  the  heat  sink.  This  is  a  very  well  researched  and  well  documented  area  in 
literature.  That  is  why  one  of  the  example  problem  compares  the  computations  of 
AMPHIB  with  some  well  established  experimental  results.  The  point  of  comparison 
is  the  velocity  field  at  the  mid-section  indicated  by  the  shaded  area.  The  velocity  field 
in  the  enclosure  is  completely  three-dimensional. 
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Figure  3.2 

The  computed  and  experimental  velocity  fields  are  compared.  The  computed  velocity 
fields  are  on  the  left.  The  experimental  ones  due  to  Arroyo  and  Saviron  (1992)  are  on 
the  right.  The  box  is  25mm  wide,  14.6mm  broad  and  12.3  mm  high.  The  fluid  is 
silicone  oil  (Pr  =  130).  From  top  to  bottom,  the  Rayleigh  numbers  are  respectively 
6319.0,  1 1 101.0,  22884.0,  44744.0  and  66604.0.  Direct  comparison  between 
experiments  and  simulations  for  three-dimensional  flows  are  extremely  rare.  This  is, 
to  the  best  of  our  knowledge,  the  first  of  its  kind. 
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Figure  3.3 

The  computed  and  experimental  velocities  are  compared  for  the  same  section  as  in 
figure  3.2.  What  is  however  being  compared  are  the  pathlines.  There  is  a  distinct 
advantage  in  comparing  pathlines,  since  they  can  be  observed  in  experiments  directly. 
On  the  other  hand,  velocity  vectors  have  to  be  deduced  from  time-lapse  photographs  in 
experiments  and  are  more  difficult  to  do.  The  experiments  are  on  the  right,  and  the 
simulations  are  on  the  left.  The  Rayleigh  numbers  starting  from  the  top  are 
6319.0,1 1 101.0  and  22884.0.  The  degree  of  agreement  between  the  two  is  striking. 
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Figure  3.4 


The  exact  dimensions  of  the  example  with  chip  and  substrate  are  shown.  This  is  the 
problem  that  was  shown  in  figure  1.2.  The  dimensions  are  specified  in  terms  of  the 
parameters  that  are  initialized  in  BLOCK  DATA  of  AMPHIB.  The  dimensions  in 
terms  of  millimeters  are  the  following:  H  =  152.,  WTH  =  203.2,  BTH  =  49.5, 
HCHIP  =  24.0,  BCHIP  =  6.0,  WCHIP  =  8.0,  BSUB  =  19.5,  YCHIP(l)  =  38.0, 
YCHIP(2)  =  76.0,  YCHIP(3)  =  1 14.0,  ZCHIP(l)  =  50.8,  ZCHIP(2)  =  101.6  and 
ZCHIP(3)  =  152.4. 
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Figure  3.5 


The  figure  at  the  top  plots  the  percentage  heat  loss  to  the  substrate  as  a  function  of  the  Rayleigh  number. 
The  Rayleigh  number  is  directly  proportional  to  the  heat  dissipated  from  the  chip.  Thus,  the  plot  shows 
that  as  the  wattage  of  the  chip  is  increased,  more  heat  as  a  percentage  is  lost  to  the  fluid  directly.  The 
figure  at  the  bottom  compares  the  average  chip  temperature  as  a  function  of  the  Rayleigh  number.  The 
linear  relationship  between  the  two  in  the  log-log  scale  is  interesting.  The  computations  carried  out 
with  AMPHIB  was  for  the  geometry  shown  in  figure  3.4  for  an  FC75  fluid  (Pr  =  30). 
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Figure  3.6 
The  figure  at  the  top  represents  flow  visualization  results  for  an  FC71  fluid  with  a 
power  dissipation  level  of  1.5W  per  chip.  The  viscosity  varies  dramatically  over  an 
order  of  magnitude.  The  enclosure  geometry  is  the  one  shown  in  figure  3.4.  The 
results  are  for  a  section  parallel  to  the  substrate  at  a  distance  of  0.5  mm  from  the  chip 
surfaces.  The  experiments  are  due  to  Joshi  et.  al  (1991).  The  computed  pathlines  for 
the  same  set  of  parameters  are  shown  at  the  bottom  figure.  The  computations  are  for 
a  fairly  fine  grid  size  of  65x38x20  and  incorporates  the  sharp  changes  in  temperature 
dependent  viscosity.  The  observed  minor  differences  in  the  flow  pattern  are  because 
the  flow  is  unsteady  and  the  exact  time  instant  has  not  been  matched.  One  observes  a 
relatively  stagnant  flow  near  the  bottom,  vertical  plumes  along  the  chips  and 
recirculating  zones  between  the  plumes  for  both  the  experiments  and  computations. 
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Figure  3.7 


The  same  case  as  seen  in  figure  3.6.  The  flow  visualization  results  are  on  the  left, 
the  computational  results  are  on  the  right.  The  results  seem  to  agree  even  in  some 
of  the  finer  details.  There  is  a  recirculating  zone  in  the  upper  left  corner.  There  is 
also  a  weak  eddy  between  the  bottom  and  the  middle  chip.  Flow  is  noticably 
stronger  in  the  upper  regions.  The  section  in  question  is  along  a  vertical  plane 
perpendicular  to  the  substrate  at  the  midplane. 
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INPUTS   AND   SUBROUTINES 


4. 1       Inputs  to  the  Program 

All  the  input  variables  are  initialized  in  BLOCK  DATA.  The  inputs  to  the  problem 
are  described  in  the  order  of  their  appearance  in  the  code. 


NIP1: 


NIP2: 
NI: 

NIM1: 


NJP1: 


NJP2: 

NJ: 

NJM1: 

NKP1 


The  number  of  control  volumes  spanning  the  x-direction  for  the  non- 
staggered  grid.  The  number  of  control  volumes  in  the  computational  domain 
will  be  two  less  since  one  control  volume  on  each  end  is  used  to  impose 
boundary  conditions. 
It  is  set  to  NIP  1  +  1. 
It  is  settoNIPl  -  1. 

It  is  set  to  NIP1  -  2.  It  is  exactly  equal  to  the  number  of  control  volumes  in 
the  computational  domain.  The  number  of  control  volumes  in  a  realistic 
simulation  depends  strongly  on  the  Rayleigh  number.  The  higher  the 
Rayleigh  number  the  greater  the  number  of  grid  points  required.  A  useful 
rule  of  the  thumb  is  that  the  average  resolution  for  natural  convection 
problems  should  at  least  0.1  (non-dimensional)  in  the  horizontal  directions 
(x  and  z)  and  at  least  0.05  in  the  vertical  direction.  Typically  for  a  0.3:1:1.3 
geometry  NIM1  is  set  to  10.  Thus,  NI,  NIP1  and  NIP2  are  respectively  11, 
12  and  13. 

The  number  of  control  volumes  spanning  the  y-direction  for  the  non- 
staggered  grid.  The  number  of  control  volumes  in  the  computational  domain 
will  be  two  less  since  one  control  volume  on  each  end  is  used  to  impose 
boundary  conditions.  NJP1  is  set  to  22  according  to  our  prescription. 
It  is  set  to  NJP1  +  1. 
It  is  settoNJPl  -  1. 
It  is  set  to  NJP1  -  2. 

The  number  of  control  volumes  spanning  the  z-direction  for  the  non- 
staggered  grid.  The  number  of  control  volumes  in  the  computational  domain 
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will  be  two  less  since  one  control  volume  on  each  end  is  used  to  impose 
boundary  conditions.  It  is  set  to  32  for  the  geometry  discussed. 

NKP2:  It  is  set  to  NKP1  +  1. 

NK:  It  is  set  to  NKP1  -  1. 

NKM1:  It  is  set  to  NKP1  -  2. 

NNMAX:        Total  number  of  variables  when  the  matrix  is  assembled.  NNMAX  is  the 
product  of  NIP1,  NJP1  and  NKP1  and  must  be  initialized  accordingly. 

NMAX:  Prescribed  number  of  time  steps  per  run.  There  is  no  convergence  criterion 

to  decide  whether  steady  state  is  reached,  since  the  algorithm  is  essentially 
integrating  forward  in  time.  NMAX  is  basically  the  only  control  on  the 
length  of  a  run.  An  NMAX  of  1000  is  typical.  In  practice  a  complete  run  is 
never  achieved  by  a  single  run  and  a  number  of  restarts  are  required. 

MAX:  The  maximum  number  of  iterations  in  the  iterative  solver,  i.e.,  even  if  the 

residual  error  criterion  is  not  met  the  number  of  iterations  does  not  exceed 
IMAX.  IMAX  has  been  set  to  10. 

ITMAX:  Maximum  number  of  iterations  in  the  pressure  loop.  Again,  if  the  mass  flux 

criterion  is  not  met,  the  number  of  iterations  cannot  exceed  ITMAX. 
ITMAX  is  set  to  5.  If  the  number  of  iterations  consistently  exceed  ITMAX, 
it  generally  means  that  the  time  step  DTIME  must  be  reduced. 

KRUN:  A  flag  to  determine  if  the  job  is  a  restart  or  a  new  run.  If  KRUN  is  set  to  0 

computations  start  from  scratch.  Otherwise,  an  input  file  containing  the  field 
variables  is  read  in  to  continue  the  computations  if  KRUN  is  set  to  1. 

NPRINT:         The  field  variables  are  written  to  a  file  every  NPRINT  time  steps.  Also, 
statistics  for  the  overall  mass  and  energy  balance  are  also  printed  every 
NPRINT  time  steps.  This  is  a  checkpointing  feature  to  ensure  that  the 
computations  are  not  wasted  if  the  program  terminates  abnormally. 
It  is  typically  set  to  100. 
SMALL:  A  variable  which  is  set  to  an  arbitrary  high  value  to  force  the  velocities  in  the 

conduction  regions  to  a  value  very  close  to  zero.  Typically  SMALL  is  set  to 
1020. 
EPS:  Minimum  prescribed  residual  in  the  iterative  solver.  In  the  subroutine  SIP 

(the  iterative  solver)  the  prescribed  square  of  the  error  norm  is  calculated 
from  EPS.  i.e.,  II  Ax  -  b  II2  <  NxEPS2,  where  N  is  the  number  of 
variables  and  the  matrix  equation  is  Ax  =  b.  A  value  of  108  with  all 
variables  in  double  precision  is  quite  typical.  The  error  norm  is  normalized 
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SORMAX: 


ISUB: 

ICHIP: 
JCHIP: 
KCHIP: 

NCHP: 


ICHOICE: 


XBR: 


YBR: 


with  the  norm  of  the  variable  unless  the  norm  is  small  (less  than  1(H)  in 
which  case  the  error  norm  is  absolute  and  not  relative.  This  is  automatically 
taken  care  of  in  the  linear  solver. 

The  maximum  permissible  absolute  sum  of  all  mass-fluxes  from  ail  control 
volumes  at  each  time  step.  This  parameter  controls  the  number  of  iterations 
in  the  pressure  loop.  In  other  words,  if  the  mass  fluxes  exceed  SORMAX 
the  outer  pressure  loop  is  traversed  again  until  the  flux  reduces  to  a  level 
below  SORMAX.  For  a  run  where  the  final  flow  field  is  expected  to  be 
steady,  the  value  of  SORMAX  can  be  fixed  quite  arbitrarily  since  at 
convergence  the  mass  fluxes  reduce  to  an  arbitrarily  small  value.  For 
unsteady  runs  (e.g.,  oscillatory  flow)  the  level  must  be  chosen  with  care. 
As  the  grid  size  is  reduced,  SORMAX  is  increased.  SORMAX  must  not 
exceed  2.0. 
The  number  of  control  volumes  spanning  the  substrate  in  the  x-direction. 

It  should  be  at  least  2. 

The  number  of  control  volumes  spanning  each  chip  in  the  x-direction. 
The  number  of  control  volumes  spanning  each  chip  in  the  y-direction. 
The  number  of  control  volumes  spanning  each  chip  in  the  z-direction. 
At  the  very  least,  ICHIP,  JCHIP  and  KCHIP  should  be  set  to  2. 
A  flag  when  set  to  0  solves  an  enclosure  problem  without  chips  or 
substrate.  When  NCHP  is  set  to  1  the  conjugate  heat  transfer  problem  is 
solved  which  now  includes  conduction  in  the  substrate  and  chip  as  well 
as  fluid  convection. 

It  is  a  parameter  which  determines  the  fluid  to  be  simulated.  If  ICHOICE  is 
set  to  0,  the  fluid  has  a  constant  viscosity.  If  ICHOICE  is  set  to  1,  the  fluid 
is  FC75  and  if  ICHOICE  is  set  to  2,  the  fluid  is  FC7 1 .  The  user  can 
simulate  any  fluid  provided  the  kinematic  viscosity  as  a  function  of  the 
temperature  is  provided.  The  changes  are  to  be  made  in  the  routines  PROP 
and  CALT. 

This  parameter  is  required  for  generating  the  grid.  It  is  the  dimension  of  the 
smallest  control  volume  in  mm  next  to  the  wall  in  the  x-direction.  The 
smallest  grid  is  always  next  to  the  wall  in  order  to  resolve  the  boundary 
layer.  Typically,  it  should  be  set  at  a  value  that  is  one  percent  of  the 
enclosure  dimension  in  that  direction  which  is  BTH  in  this  case. 
The  smallest  grid  size  in  the  y-direction  in  mm.  Again,  one-hundredth 
of  H  is  suggested.  However,  for  a  general  case,  sensitivity  to  these 
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ZBR: 

H: 

WTH: 

BTH: 

HCHIP: 

WCHIP: 

BCHIP: 

BSUB: 

ALC: 

ALS: 

THOT: 

TCOOL: 

TAVG: 


RHC: 

RHS: 

QQQ: 
DTIME: 


PR: 


RA: 


parameters  should  be  tested. 

The  smallest  grid  size  in  the  z-direction  in  mm.  A  similar  prescription  holds. 

The  height  of  the  enclosure  in  mm.  This  is  the  y-direction. 

The  width  of  the  enclosure  in  mm.  This  is  the  z-direction. 

The  length  of  the  enclosure  in  mm.  This  is  the  x-direction. 

The  dimension  of  an  individual  chip  in  the  y-direction  in  mm. 

The  dimension  of  an  individual  chip  in  the  z-direction  in  mm. 

The  dimension  of  an  individual  chip  in  the  x-direction  in  mm. 

The  thickness  of  the  substrate  in  mm  (x-direction). 

Thermal  conductivity  ratio  between  chip  and  fluid. 

Thermal  conductivity  ratio  between  substrate  and  fluid. 

The  non-dimensional  temperature  at  the  bottom  wall. 

The  non-dimensional  temperature  at  the  top  wall. 

The  reference  temperature  in  °C.  This  is  required  to  evaluate  the  viscosity  as 

a  function  of  the  temperature.  In  the  case  when  the  top  and  bottom 

temperatures  are  equal,  THOT  and  TCOOL  is  set  to  zero,  the  TAVG  is 

set  equal  to  the  temperature  of  the  top  and  bottom  walls. 

The  heat  capacity  (product  of  density  and  specific  heat)  ratio  between  chip 

and  fluid. 

The  heat  capacity  (product  of  density  and  specific  heat)  ratio  between 

substrate  and  fluid. 

This  is  the  heat  dissipated  per  chip  in  Watts. 

Non-dimensional  time  step.  The  explicit  time  step  for  an  upwind  scheme 

is  calculated  by  the  routine  TSTEP.  Since  the  scheme  used  for  the  program 

is  QUICK,  it  is  recommended  that  the  time  step  is  about  a  third  of  the  time 

step  calculated. 

Prandtl  number  of  the  fluid.  This  is  a  non-dimensional  property  of  the  fluid 

defined  as  — ,  where  v  is  the  kinematic  viscosity  and  a  is  the  thermal 
a 

diffusivity. 

The  product  of  the  Prandtl  number  and  the  Rayleigh  number  which  is 

eSqH5 

^-* — ,  where  g  is  the  acceleration  due  to  gravity;  p  is  the  coefficient  of 

a2Kc 

volume  expansion;  L  is  the  length  scale  (height  of  the  enclosure  in  this 
case);  and  Kc  is  the  thermal  conductivity  of  the  chip;  q  is  the  volumetric  rate 
of  heat  generation  for  each  chip.  This  must  be  specified  only  if  the  it  is  a 


36 


pure  enclosure  problem  (when  NCHP  =  0).  Otherwise,  the  product  of  the 
Rayleigh  and  Prandtl  numbers  are  automatically  calculated  using  QQQ. 

IUNFRM:        This  parameter  is  used  in  the  grid.  When  set  to  0,  a  non-uniform  grid  is 
generated.  Otherwise,  a  uniform  grid  is  created. 

QUICK:  It  is  a  number  between  0  and  1 .  If  QUICK  is  set  to  0.  the  upwind  scheme 

is  used.  If  QUICK  is  set  to  1.  the  quick  scheme  is  used.  Any  number 
between  the  two  is  a  weighted  average  of  the  two  schemes.  It  is 
recommended  that  for  initial  runs  especially  for  high  Rayleigh  numbers 
(greater  than  107)  QUICK  be  set  to  0.  After  a  converged  solution  is 
got  QUICK  should  be  set  to  1.  However,  if  there  are  convergence  problems 
QUICK  should  be  set  to  zero  again.  A  failure  to  converge  to  a  steady  state 
(when  experiments  clearly  indicate  steady-state)  implies  that  the  number 
of  grid  points  are  insufficient.  Therefore,  NI  etc  must  be  increased. 

NJCHIP:  The  number  of  chips  in  the  vertical  y-direction. 

NKCHIP:        The  number  of  chips  in  the  vertical  z-direction. 

YCHIP:  A  one-dimensional  array  of  length  NJCHIP.  These  are  the  locations  of  the 

chip  centers  in  the  y-direction. 

ZCHIP:  A  one-dimensional  array  of  length  NKCHIP.  These  are  the  locations  of  the 

chip  centers  in  the  z-direction. 

The  schematic  of  a  typical  geometry  that  is  solved  is  shown  in  figure  1.1.  The 
substrate  and  the  chip  protrusions  are  shaded.  In  terms  of  the  variables  in  the  program,  the 
geometry  is  marked  in  figure  3.4  as  an  example  for  a  3  by  3  three  array.  In  the  most  general 
situation,  the  protrusions  are  located  in  a  non-uniformly  spaced  NJCHIP  by  NKCHIP 
array.  All  protrusions  have  the  same  dimensions,  which  is  HCHIPx  WCHIPx  BCHIP. 


4.2       Subroutines 

PROGRAM  AMPHIB: 

This  is  the  main  program.  It  consists  of  five  subroutine  calls  in  sequence.  The  first 
subroutine  call  is  for  the  subroutine  OPENF. 
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SUBROUTINE  OPENF 

This  subroutine  opens  the  files  that  handle  the  output  and  input  for  the  code.  The  files  are 

the  following: 

UNIT  8:  This  is  the  input  data  for  the  field  variables  (tod(i,j,k),  uod(i,j,k),vod(i,j,k), 

wod(i,j,k)  and  pod(ijjc))  i.e.,  the  initial  temperature,  velocities  and 

pressure.  The  field  variables  are  read  in  only  if  KRUN  is  set  to  1. 
UNIT  10:         This  is  an  output  file  of  the  u,v  and  w  velocities  at  a  specific  location  for 

all  the  time  steps.  This  was  used  to  gauge  the  dynamical  behavior  of  the 

system  i.e.,  whether  it  is  steady-state,  oscillatory  etc. 
UNIT  1 1 :         This  is  output  file  for  the  field  variables  updated  at  specific  instances 

which  is  determined  by  NPRINT.  This  is  a  checkpointing  procedure. 
UNIT  12:         This  is  the  output  file  for  the  mean  field  variables  over  the  entire  run, 

i.e.,  over  NMAX  time  steps. 
UNIT  13:         This  is  the  file  where  the  grid  for  the  given  run  is  stored. 

SUBROUTINE  GRID: 

This  subroutine  generates  the  three-dimensional  grid  for  the  problem.  The 
subroutine  non-dimensionalizes  the  dimensions  first  in  terms  of  the  height  of  the  cavity,  h. 
It  is  therefore  sufficient  to  input  the  data  in  consistent  units  (say  inches)  rather  than  in 
millimeters.  The  subroutine  generates  a  non-uniform  grid  x  of  the  following  type 

^  _  Tanh  (HKr)  (4  1} 

K. 

This  is  known  as  the  Robert's  transformation  (Robert,  1970).  The  grid  distribution  ensures 
that  the  boundary  layers  are  resolved  near  the  wall.  At  the  same  time,  the  grid  is  spaced  as 
uniformly  as  possible  away  from  the  wall.  There  are  two  grid  parameters  (H  and  K)  that 
need  to  be  determined.  The  variable  r  represents  the  uniform  grid.  Equation  4.1  thus  maps 
a  uniform  grid  (r)  to  a  non-uniform  one  (x). 

The  grid  distribution  generated  for  the  NJCHIP  by  NKCHIP  enclosure  is  a  little 
more  complicated.  The  non-uniform  hyperbolic  tangent  distribution  of  the  grid  is  confined 
to  a  region  near  the  walls.  Otherwise,  the  grid  is  made  as  uniform  as  possible,  away  from 
the  wall.  The  same  thing  applies  to  the  y  and  z  directions.  As  an  example,  consider  a  3  by  3 
array  that  was  portrayed  in  figure  1.2.  Two  perpendicular  sections  of  the  non-staggered 
control  volumes  for  a  20x38x65  grid  is  shown  in  figure  4.1.  The  x-y  plane  is  on  the  left 
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and  the  z-y  plane  is  shown  on  the  right.  The  grid  point  is  located  at  the  centroid  of  the 
control  volume  for  the  non- staggered  grid. 

SUBROUTINE  GRID  1 

This  subroutine  is  called  by  GRID  to  calculate  the  grid  parameters  for  the  non- 
uniform grid  (H,  K).  The  routine  solves  transcendental  equations  by  the  bisection  and  the 
Newton-Raphson  method.  Note  that  the  users  could  supply  their  own  grids  in  which  case 
GRID  and  GRID1  are  not  required.  The  whole  program  can  only  handle  grids  that  are 
rectangular.  The  non-uniformity  is  restricted  to  their  respective  directions  only. 

SUBROUTINE  PROP 

In  this  routine  the  conductivity  ratios  and  the  heat  capacity  ratios  are  assigned  for 
the  chip,  substrate  and  fluid.  The  heat  capacity  is  the  product  of  the  specific  heat  and 
density.  The  parameter  ICHOICE  determines  the  fluid  chosen.  In  the  present  program,  the 
choice  is  limited  to  FC71  and  FC75.  Using  the  properties,  the  product  of  the  Rayleigh  and 
Prandtl  numbers  (sometimes  referred  to  as  the  Boussinesq  number)  is  calculated  using 
QQQ,  the  heat  dissipated  per  chip.  The  user  will  need  to  modify  this  subroutine  accordingly 
for  a  different  set  of  fluids. 

SUBROUTINE  INITIO 

In  this  subroutine  all  variables  are  initialized.  If  the  computations  are  started  from 
scratch,  the  field  variables  are  all  initialized  to  zero.  For  a  restart  job  (KRUN  =  1),  the  field 
variables  are  read  in  from  UNIT  8.  If  KRUN  is  set  to  0,  and  NCHP  set  to  zero  as  well,  a 
problem  of  pure  natural  convection  is  solved. 

SUBROUTINE  PLOOP 

This  subroutine  incorporates  the  essence  of  the  SIMPLEX  algorithm.  It  also 
includes  an  error  control  routine  (Liu,  1979).  It  calls  the  following  subroutines: 

SUBROUTINE  CALT 

This  subroutine  calculates  the  coefficients  of  the  algebraic  equation  for  solving  the 
finite-difference  equation  for  temperature  which  is  of  the  following  form: 

ApTp=  AETE  +  AWTW  +  ANTN  +  ASTS  +  AETE  +  ABTB 

+  AFTF  +    S  (4.2) 
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The  explicit  expressions  for  the  coefficients  (i.e.,  Apetc)  as  well  as  the  derivation  of  the 
equations  from  scratch  are  given  in  appendix  B  for  the  interested  reader. 

The  coefficients  must  be  modified  for  accommodating  the  boundary  conditions. 
This  is  also  accomplished  in  the  subroutine.  The  computational  details  are  available  in 

appendix  B. 

The  interfacial  thermal  conductivities  are  calculated  in  this  routine.  It  uses  the 
harmonic  mean  formulation  as  suggested  by  Patankar  (1980).  Essentially  the  conductivities 
are  calculated  in  the  following  manner: 

Axi  +  Axi+1  r4~ 

Ke  =  Ki  Ki+1  (- )  (4.3) 

AXiKi+i  +  Axi+1Kj 

Ke  is  the  thermal  conductivity  of  the  east  face.  In  this  manner,  sharp  changes  in 
conductivities  can  be  handled  accurately.  The  formulation  is  necessary  since  properties  vary 
dramatically  in  a  conjugate  heat  transfer  problem  such  as  this.  For  instance,  in  a  problem 
with  aluminum  protrusions,  plexiglas  substrate  and  FC75  liquid  the  conductivity  ratio 
between  aluminum  and  FC75  is  almost  2800.  For  the  heat  capacity  ratios,  the  interfacial 
heat  capacities  are  calculated  by  linear  interpolation. 

SUBROUTINE  CALU 

The  finite  difference  analog  of  the  x-momentum  equation  is  assembled  in  this 
subroutine.  The  equations  are  linearized  and  the  unknown  to  be  solved  is  the  U-velocity  at 
every  grid  point.  The  equation  is  of  the  following  form: 

ApUp=  AEUE  +  AWUW  +  ANUN  +  ASUS  +  AEUE  +  ABUB 

+  AFUF  +    S  (4.4) 

The  boundary  conditions  are  also  incorporated  by  suitable  modifications  of  the  coefficients. 
All  details  are  provided  for  in  appendix  B.  The  equation  for  solving  the  coupling  constants 
(to  be  discussed  later)  of  the  U-velocity  for  each  grid  point  is  also  calculated.  The  coupling 
constants  will  be  used  later  in  the  equation  for  the  pressure  correction. 

Similarly,  the  subroutines  CALV  and  CALW  are  used  to  assemble  the  equations 
and  apply  boundary  conditions  for  the  V  and  W  velocities  respectively.  The  equation  for 
solving  the  coupling  constants  are  also  assembled  in  these  two  subroutines. 

In  all  the  velocity  routines,  the  interfacial  kinematic  viscosity  for  the  fluid  region  is 
calculated  by  a  harmonic  mean  formulation. 
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SUBROUTINE  CALP 

The  pressure  correction  equation  is  derived  from  the  continuity  equation  and 
enforces  mass  balance  by  correcting  the  pressure  and  velocities.  The  fundamental 
relationships  are  the  following: 

U*e  =  U'-due(P'E-P'p)  (4-5) 

V*n  =  V'-dvn(P'N-P'p)  (4-6> 

W*f  =  W  -  dwf  (P'F  -  P'p)  (4-7) 

The  starred  quantities  are  corrected  velocities  that  satisfy  the  continuity  equations.  The 
primed  quantities  for  pressure  are  the  pressure  correction.  The  coupling  constants  due,  dvn, 
and  dwf  were  determined  respectively  in  CALU,  CALV  and  CALW.  A  pressure  correction 
equation  similar  in  form  to  equation  (6)  is  set  up  in  CALP. 

SUBROUTINE  NU 

This  subroutine  calculates  the  average  Nusselt  numbers  at  the  cold  and  hot  walls 
and  the  generation  of  internal  energy.  It  does  an  overall  heat  balance.  The  mass  balance  and 
the  heat  balance  statistics  are  printed  in  this  subroutine.  NU  is  called  by  PLOOP  every 
NPRINT  time  steps. 

SUBROUTINE  SIP 

This  subroutine  is  called  by  CALU,  CALV  and  CALW  twice  every  time  step  for 
solving  the  generic  linear  equations  that  have  been  assembled  (equation  2.2).  The  first  call 
is  to  solve  the  velocity  equations  and  the  second  call  is  made  to  solve  for  the  coupling 
constants  as  dictated  by  the  SIMPLEX  algorithm.  SIP  is  called  by  CALT  and  CALP  once 
to  solve  the  temperature  and  pressure  correction  equations  respectively.  SIP  is  an  iterative 
solver.  It  is  controlled  by  two  parameters  EPS  and  IMAX.  EPS  determines  the  relative 
error  norm.  EPS  is  generally  set  to  108  for  double  precision  calculations.  Double  precision 
calculations  are  recommended  and  are  in  fact  the  default  (using  the  IMPLICIT  statement  in 
every  subroutine). 

IMAX  is  the  maximum  number  of  iterations  allowed  in  the  iterative  solver.  IMAX 
is  typically  set  to  a  value  of  10.  The  subroutine  calls  the  subroutines  RES  and  XL.  RES 
calculates  the  error  vector  i.e.,  Ax  -  b.  XL  performs  the  forward  and  back  substitution  for 
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the  incomplete  LU  decomposition  which  is  an  integral  part  of  the  solver.  The  solver  is 
known  as  Strongly  Implicit  Procedure  and  is  credited  to  Stone  (1968).  The  reader  can 
substitute  any  solver  in  its  place.  SIP  was  found  to  be  the  most  cost-effective  for  natural 
convection  problems  that  were  studied.  Note  that  iterative  solvers  are  generally  more 
efficient  than  direct  solvers  for  three-dimensional  problems  since  the  matrices  that  are 
assembled  are  generally  speaking  ill-conditioned.  The  condition  number  (ratio  of  the  largest 
to  the  smallest  eigen  value  of  the  matrix)  increases  for  finer  grids.  If  the  smallest  eigen 
value  of  the  matrix  is  exacdy  zero,  then  the  matrix  is  singular. 

The  input  (arguments)  to  the  subroutine  are  1ST,  JST,  KST,  ISPJSP  and  KSP  and 
the  variable  to  be  solved  (the  field  variable).  1ST,  JST  and  KST  are  the  starting  indices  of 
the  computational  domain  for  the  x,  y  and  z  direction.  ISP,  JSP  and  KSP  are  the  last 
indices  for  the  respective  directions. 

SUBROUTINE  TSTEP 

This  subroutine  calculates  the  maximum  allowable  time  step  for  the  difference 
equations  of  the  temperature  and  the  velocities  if  the  scheme  were  explicit.  This 
information,  i.e.,  the  minimum  time  step  is  printed  every  NPRINT  time  steps.  This  routine 
is  called  by  PLOOP.  This  provides  a  useful  guideline  for  the  time  step  that  is  specified  by 
the  variable  DTIME.  The  explicit  time  step  is  calculated  using  the  CFL  criteria  and  for  an 
upwind  interpolation  scheme.  The  criteria  for  the  QUICK  scheme  will  be  more  stringent.  It 
is  therefore  recommended  that  the  time  step  be  assigned  a  value  that  is  between  a  third  or 
half  the  computed  time  step.  In  this  way,  numerically  induced  oscillations  can  be  avoided. 

SUBROUTINE  CHIPTEMP 

This  subroutine  calculates  the  average  temperature  rise  of  the  chips.  The  temperature  is 
calculated  by  volume  average  of  temperatures  of  the  control  volumes  that  are  contained  in 
the  chip.  The  average  temperatures  are  then  suitably  scaled  to  give  the  temperature  rise  in 
degree  Celsius.  The  relevant  information  is  then  printed.  The  subroutine  is  called  by 
PLOOP  every  NPRINT  time  steps. 
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Figure  4.1 


The  figure  on  the  left  and  right  are  two  perpendicular  views  of  the  grid 
generated  for  a  fine  mesh  (20x38x65).  The  control  volume  sections  of 
the  non-staggered  grid  points  are  what  is  shown.  The  grid  generated 
corresponds  to  a  3x3  array  shown  in  figure  3.4.  Some  of  the  results 
produced  for  the  above  grid  are  shown  in  figures  3.6  and  3.7. 
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Figure  4.2 


This  is  a  schematic  flow  chart  of  the  subroutines  in  the  code  AMPHIB. 
The  main  program  has  five  subroutine  calls;  OPENF,  GRID,PROP, 
INITIO  and  PLOOP.  The  subroutine  PLOOP  is  blown  up  to  reveal 
more  details. 
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APPENDIX  A 


5.1        Boundary  Conditions 
5.1.1    The  Defaults 

The  temperature  boundary  conditions  have  been  set  up  in  the  program  as  follows: 
Top  wall:         Isothermal  (perpendicular  to  y  axis) 
Bottom  wall:     Isothermal  (perpendicular  to  y  axis) 
Left  wall:  Adiabatic  (perpendicular  to  x  axis) 

Right  wall:  Adiabatic  (perpendicular  to  x  axis) 
Back  wall:  Adiabatic  (perpendicular  to  z  axis) 
Front  wall:  Adiabatic  (perpendicular  to  z  axis) 
The  velocity  boundary  conditions  are  all  non-slip. 

The  temperature  boundary  conditions  are  imposed  in  the  subroutine  CALT.  The 
scope  of  the  program  can  be  expanded  by  imposing  different  combinations  of  boundary 
conditions.  The  three  commonly  used  boundary  conditions  are:  (1)  Adiabatic  (2) 
Isothermal,  (3)  Convective  or  mixed  and  (4)  specified  heat  flux.  Each  of  these  will  now  be 
described  in  detail.  The  velocity  boundary  conditions  imposed  in  the  program  are  all  non- 
slip. 

5.1.2    Adiabatic  Boundary  Conditions 

This  corresponds  to  the  condition  of  having  the  heat  flux  set  to  zero  at  the  wall.  For 
example,  V-  (at  the  left  wall)  =  0.  The  numerical  boundary  condition  corresponding  to 

this  for  the  control  volume  next  to  the  left  wall  (i  =  2)  is  Tw  =  TP.  Substituting  this 
condition  in  the  algebraic  equation  (equation  2.2),  the  coefficients  must  be  modified  in  the 
following  manner: 

Ap  =  Ap  -  Aw  and  Aw  =  0. 
For  the  control  volume  that  is  one  away  from  the  wall  (i  =  3),  the  numerical  boundary 
condition  is  Tww  =  Tw-  This  corresponds  to  S  =  Ty/Aww  in  the  program. 
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5.1.3  Tsothermal  Boundary  Conditions 

This  corresponds  to  the  condition  where  the  temperature  at  the  wall  is  specified.  For 
example,  T(at  the  bottom  wall)  =  THOT  .  The  numerical  boundary  condition  corresponding 
to  this  for  the  control  volume  next  to  the  bottom  wall  (j  =  2)  is  Ts  +  Tp  =  2  THOT. 
Substituting  this  condition  in  the  algebraic  equation  (equation  2.2),  the  coefficients  must  be 
modified  in  the  following  manner: 

AP  =  AP  +  As  ,  S  =  2  THOT  As  and  As  =  0. 
For  the  control  volume  that  is  one  away  from  the  wall  (j  =  3),  the  numerical  boundary 
condition  is  TSS  =  2  THOT  -  Ts.  This  corresponds  to  S  =  (2  THOT  -  Ts)  AWw  in  the 
program. 

5.1.4  Convective  Boundary  Condition 

This  corresponds  to  the  case  when  the  temperature  at  the  wall  satisfies  a  boundary 
condition  of  the  type  |£  +  (HCON)  T  =  0,  for  example  at  the  back  wall.  The  numerical 
boundary  condition  corresponding  to  this  for  the  control  volume  next  to  the  bottom  wall 

(k=2)  is  HCON(Tp  i  Tfi  )  +  Tp  "  Tb  =  0.  Substituting  this  condition  in  the  algebraic 
2  Az 

equation  (equation  2.2),  the  coefficients  must  be  modified  in  the  following  manner: 

1        H 

Az       l 
AP  =  AP  -  CONST  AB  ,  where  CONST  =  — ^ 

Az      2 
For  the  control  volume  that  is  one  away  from  the  wall  (k  =  3),  the  numerical  boundary 
condition  is  TBb  =  CONST  TB.  This  corresponds  to  S  =  (CONST  )TB  ABb  in  the  program. 

5.1.4    Heat  Flux  Boundary  Condition 

This  corresponds  to  the  condition  of  having  the  heat  flux  specified  at  the  wall.  For 
example,  V-  (at  the  left  wall)  =  C.  The  numerical  boundary  condition  corresponding  to 

this  for  the  control  volume  next  to  the  left  wall  (i  =  2)  is  Tw  =  TP  -  Const,  where  Const  = 
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CAx.  Substituting  this  condition  in  the  algebraic  equation  (eqn  5),  the  coefficients  must  be 
modified  in  the  following  manner: 

Ap  =  Ap  -  Aw  ,  S  =  -  Aw  Const  and  Aw  =  0. 
For  the  control  volume  that  is  one  away  from  the  wall  (i  =  3),  the  numerical  boundary 
condition  is  TWw  =  Tw  -  Const.  This  corresponds  to  S  =  (Tw  -   Const )  Aww  in  the 
program.  Note  that  the  adiabatic  condition  is  a  special  case,  when  C  =  0. 

5.2       Velocity  Boundary  Conditions 

All  velocities  i.e.,  the  U,  V  and  W  velocities  are  set  to  zero  (non-slip  condition)  for 
all  six  walls  in  the  program.  The  most  commonly  used  boundary  conditions  and  their 
implementation  will  now  be  described. 

5.2.1  Non-Slip  Conditions 

The  numerical  implementation  is  complicated  by  the  staggered  nature  of  the  grid. 
Consider  the  left  wall  as  an  example.  For  the  u-velocity  the  grid  point  coincides  with  the 
left  wall  due  to  stagger).  Hence,  in  equation  2.2  we  simply  set  Aw  =  0  (since  Uw  =  0). 
For  the  V  and  W  velocities  for  the  control  volume  next  to  the  wall  (i  =  2  for  V  and  W,  and  i 
=  3  for  U)  we  have  VP  =  -  Vw  and  WP  =  -  Ww.  Hence  AP  =  AP  +  Aw  and  Aw  =  0  in  each 
of  the  cases. 

For  the  control  volume  one  away  from  the  wall  (i  =  3  for  V  and  W,  and  i  =  4  for  U) 
we  have  Uww  as  identically  zero.  Hence  Aww  =  0  for  the  U-velocity  grid.  For  the  V  and  W 
velocities  we  have  Vw  =  -  Vww  and  Ww  =  -  Www-  Hence,  S  =  -  AwwVw  and  Aww  =  0  for 
the  V  velocity  equation  and  S  =  -  AwW\v  and  Aww  =  0  for  the  W  velocity  equation. 

5.2.2  Slip  Wall  Conditions 

Again  consider  the  left  wall  as  an  example.  These  boundary  conditions  are 
employed  to  generally  impose  a  symmetry.  For  instance,  if  the  flow  consists  of  two 
symmetric  vortices,  a  slip  wall  can  be  imposed  on  one  of  them  and  calculations  can  be 
carried  out  for  only  half  the  computational  domain.  The  slip  wall  is  the  numerical 
implementation  of  zero  shear  stress  at  the  wall.  For  the  left  wall  it  would  mean  that  -^  =  0 


47 


and  <W   _  Q  AlsQ  u  =  q  Tne  numerical  boundary  conditions  for  the  U  velocity  is  the 
dx 

same  as  the  non-slip  conditions  at  the  wall  i.e.  Aw  =  0  (i  =  3)  and  AWw  =  0  (i  =  4).  The 
boundary  conditions  for  the  V  and  W  are  identical  to  the  adiabatic  condition  for  the 
temperature  i.e.,  AP  =  AP  -  Aw  and  Aw  =  0  (i  =  2)  and  S  =  TwAww  and  Aww  =  0  (i  =  3). 

5.2.3    Velocity  Prescribed 

This  is  a  very  important  boundary  condition.  Using  this  one  can  solve  problems  in 
mixed  or  forced  convection  as  in  figure  1.3.  As  an  example,  let  the  left  wall  have  an  inlet 
with  a  U-velocity  of  U0  (consistently  non-dimensionalized).  Also  the  V  and  W  velocity  are 
zero.  For  that  portion  of  the  wall,  the  boundary  conditions  for  the  U  velocity  are  S  =  UoAw 
and  Aw  =  0  (i  =  3)  and  S  =  UoAww  and  Aww  =  0  (i  =  4).  The  V  and  W  velocities  have  the 
usual  non-slip  boundary  condition. 

For  the  control  volume  next  to  the  wall  (i  =  2  for  V  and  W,  and  i  =  3  for  U),  we 
have  VP  =  -  Vw  and  WP  =  -  Ww.  Hence  AP  =  AP  +  Aw  and  Aw  =  0  in  each  of  the  cases. 
For  the  control  volume  one  away  from  the  wall  (i  =  3  for  V  and  W),  we  have  Vw  =  -  Vww 
and  Ww  =  -  Www-  Hence,  S  =  -  AwwVw  and  Aww  =  0  for  both. 

If  there  is  an  inlet  there  must  be  an  outlet.  For  the  outlet  the  velocities  cannot  be 
prescribed.  In  fact,  only  approximate  boundary  conditions  can  be  used.  This  is  the  so 
called  'natural  boundary  condition'.  For  instance,  if  the  outlet  is  at  a  portion  of  the  east  wall 
we  numerically  implement  ^  =  0,    -^  =  0  and   -^-  =  0  at  the  east  wall.  The  numerical 

implementation  of  these  conditions  have  already  been  discussed  and  are  quite 
straightforward.  Basically,  this  means  that  we  are  not  sure  what  the  boundary  conditions 
are  but  assume  that  they  don't  change  much  near  the  outlet.  Overall  mass  balance  is 
automatically  satisfied. 
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APPENDIX   B 


6. 1       Difference  Equations 

In  this  section  the  finite  difference  equations  are  presented  for  the  equations  of  the 
velocities  and  temperatures.  The  details  of  the  derivation  are  skipped.  Only  the  final  form  of 
the  equations  are  stated.  For  additional  details  of  the  control  volume  method  the  reader  is 
referred  to  Patankar  (1980). 

6.1.1    Energy  Equation  and  linear  equation  for  Temperature 

The  QUICK  scheme  for  estimating  temperature  fluxes  at  the  east  and  west  faces  are 
the  following: 

n  ^  r    {  Axi+1TP  + AxjTe 

Oe  1  e     =     <Je  \  2AXe 

.  JGgj  -  Ge    w  AXjAXj-n  TpAXee  TEEAxe  y     , 

+      K  8  M  AXgAXee  M  Axe  +  Axee  Axe  +  Axee 

lGel  +  Ge  w  AxjAxi+i  TEAxw  TwAxe        _T)  (61) 

"     ^  8  M  AxeAxw  MAxw  +  Axe  Axe  +  Axw 

n   t  r     (  Axj.iTp  +  AxjTw 

Uwlw    -    ^Jwl  2Axw  ' 

1GWI  -  Gw  w  AxjAxj.i          TEAxw  TwAxe        _T 

+    ^           8          M  AxeAxw  M  Axe  +  Axw  Axe  +  Axw          PJ 

IGWI  +  GW        AxjAxj.i  TwwAxw  TpAx^ _  T    ^    (6  2\ 

'    \           8           M  AxwwAxw  M  Axw  +  Axww  Axw  +  Axww 

The  coefficients  of  the  energy  equation  are  the  following: 

GeAxj  AyjAzk  IGel  +  Ge  AxjAxi+1 

Ae  =    '  ~2Ax7    +       Axe       +     {  8  M  Axe(Axw  +  Axe) 

,    IGel  -Ge    wAXjAxi+i  1GWI  -  Gw  AXjAXj.! (63) 

+      < 8 )(  AxeAxee)      +       (_      8  }  Axe(Axw  +  Axe)  ^ 
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GwAxi  AyjAzk  , 1GWI ■  +  Gw  >.  ,  AxjAxj.i 

Aw=      "2Ax7    +       Axw  (  8  MAxwAxww; 

.  1GWI  -  Gw  x        AxjAxj.!  rIGel  +  Ge  , AxjAxi+1 (64) 

+    C  g  )  Axw(Axw+  Axe)         ^8  '  Axw(Axw  +  Axe) 

GnAyi  AxjAzk  IGnl  +  Gn       AyiAyi+i 

An  =    "  ~2Ay7  ~Ay7~  (  8  M  Ayn(Ays  +  Ayn) 

IGnl  -  Gn       AyjAy^  (  IG5'- Gs  )  A     %&V      ,  (6.5) 

+    <  8  HAynAynn;  ^         8         ;  Ayn(Ays  +  Ayn) 

GsAyi  AxjAzk  IGSI  +  Gs       AyjAyj.u 

A*  =      "2Ay7    +     "^yT  (         8  HAysAyss; 

IG,|  -  Gs  AyiAyH (!Gnl_t^)^^Ayi±L_^  (6.6) 

+    (         8         }  Ays(Ays+  Ayn)     +   ^     8  ;  Ays(Ays  +  Ayn) 

GfAzk  AxjAyj  IGfl  +  Gf  AzkAzk+i 

Af  =    "  "2Az7    +       Azf  K         8  ]  AzKAzb  +  Azf) 

,  IGfl  -  Gf  w  AzkAzk+1  IGbl  -  Gb  , AzkAzk.! (6  ?) 

+    ( 8 )(   AzfAzff    }    +     (         8  )  Azf(Azb  +  Azf)  V°'/J 

GbAzk  AxjAy,  IGbl  +  Gb  ,  AzkAzk.! 

Ab  "        2Azb     +       Azb  K         8  ;  AzbAzbb 

,  IGbj  -  Gh  ,        AzkAzk.i  IGfl  +  Gf,    AzkAzk+1 6  g) 

+    < 8  }  Azb(Azb  +  Azf)    +   K       8  ;  Azb(Azb  +  Azf)  ^'  ; 

_        AxjAyjAzk  rr  n-i        ,  IGel  -  Ge  x  AxjAxj+i  T 

S=  At         lp      "    (         8  }  Axee(Axee  +  Axe)    lEE 

IGWI  +  Gw  AxjAxj.! 

'    K  8  ;  Axww(Axww  +  Axw)     w 

,  IGnl  -  Gn  ,  AyiAyi+1 

"    (        8         }  Aynn(Aynn  +  Ayn)    1nn 
,  IGSI  +  Gs  AyjAyj.i  y 

"     (  8  J  Ayss(Ayss+ Ays)    lss 

,  IGfl  -  Gf  .      _AzkAzk±i__  T 
'    ^         8         J  Azff(Azff+ Azf)    lFF 
IGbl  +  Gb  ,  AzkAzk4 

•    (  8  }  Azbb(Azbb  +  Azb)    1bb  (6-9) 


AP=      AE  +    Aw  +    AN  +     As+     AF  +    Ab 

IGel  -  Ge    N  AXjAXj+i 

*  8         ;  Axee(Axee  +  Axe) 

IGWI  +  Gw  v  AxjAxj.i 

'    {  8  ;  Axww(AxWw  +  Axw) 

,   IGnl  -  Gn  AyiAyi+1 

v  8         ;  Aynn(Aynn  +  Ayn) 

(  IG,I  +  Gs  AyiAyM 

^  8  >  Ayss(Ayss+ Ays) 

f  IGfl  -  Gf  AzkAzk+i 

1         8         ;  Azff(Azff+ Azf) 

.  IGbl  +  Gb  v  AzkAzk.i  AxiAyiAzk  ffi  im 

"     (  8  >  Azbb(Azbb  +  Azb)  At  v  •     } 
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The  coefficient  Ap  etc  have  their  usual  meaning. 

6.1.2    x-Momentum  Equation  and  Linear  Equation  for  u 

For  the  x-momentum  equation,  the  velocities  must  be  interpolated  differendy  due  to 
the  staggered  nature  of  the  grid.  The  QUICK  scheme  results  in  the  following  for  the  east 
and  west  faces  of  the  control  volume: 

Up  +  uE    x 
CjeUe      =      <Je  ( o  > 


2 
C 
16        M  Axe  }  v  Axi+i  Axi 

2 

W   ~ }  v  Ax~  '  v    Axi  Axi.i 


.   IGel  -  Ge  v  ,  Axj2      uEE  -  Ue        Ue-  up  x 
+    (         Ta        m  Av  _>  v   Ay:,,  aX;      ) 


IGel  +  Ge        AXj2       uE^_Up_      jlpJ^Uv^,  v  (611) 

-    v         i6         M  Ax    M  Ax:  ' 

Up  +  UW      v 

Owuw    =    uw  v       2         ' 


.  IGWI  -Gw  x  ,  Ax;,!2      uE-  up  up-  uw   x 

+    (  "       16         M  Axw   M  Axi  AxM     ; 

2 

16~    ~  '  k  Axww  '  K  Axm  Axj.2 


,  IGWI  +  Gw  x  /Axj.i      ,uP-  uw_      uw  uWw  v  (6  l2\ 


The  coefficients  of  the  discretized  x-momentum  equations  are  the  following: 

Ge  pAVjAZk       ,  IGel  +  Ge    ,    AXj 

Ae=    "  T    +     ^   Ax~  (         16         }Ax, 


'W 

IGel  -  Ge  ,   Axj  ,  IGWI  -  Gw  x    Ax..!2 


IUel  -  ue    .     ttXj  .vjw'-^w    x     ^i-l  (6  13) 

+    ^         8         ;Axi+i  *■         16         ;AxwAXl 

Gw          r^AvjAzk  IGW1  +  GW  s  Axj.i 

Aw  =      X    +     *  AxTT  +     (          8           }  Ax,2 

,  IGWI  -  Gw  .  Axj.i  IGel  +  Ge  ,      Ax;2                                   f614. 

+    (         16         }Axw  *       16         ^AxwAxm                                l°-     ) 

a  GnAyi  AxwAzk  IGnl  +  Gn  AyiAyi+i 

An=    "iKyt  ^^A^T  ("       8           )(Ayn(Ays+Ayn) 

IGnl  -  Gn       AyiAyi+1  IGSI  -  Gs     AyjAy^ (6  15) 

+    ( 8  )(AynAyt:)  +     (         8         }  Ayn(Ays  +  Ayn)       &A:» 

a  GsAyi         rAxwAzk  IGSI  +  Gs        Ay^u 

As  =       2Ays     +    "  2Ays       +      C        8  M  AysAy    ; 


IG,I  -  Gs  AyjAyj.!  IGnl  +  Gn     Ay,Ayi+1 6 

+    (         8         }  Ays(Ays+  Ayn)    +e     8  }  Ays(Ays  +  Ayn)    {(3A0) 
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GfAzk  _  AxwAyj  IGfl  +  Gf  AzkAzk+i 

=    "  "^Az7    +  2Azf      +     ^         8  ;  Azf(Azb  +  Azf) 

IGfl  -  Gf  NAzkAzk+i  IGbl  -  Gb  ,        AzkAzk.i 


2Azf     T     il   2Azf 
IGfl  -  Gf    AzkAzk+i 
+    (-"8 }  AzfAzff       +     (  8  ^zKAzb+Azf) 


GbAzk  AxwAyj  IGbl  +  Gb     AzkAzk.! 

*  ~        2Azb     +  2Azb  l         8  ;  AzbAzbb 

,  IGbl  -  Gb  ,        AzkAzk.!  IGfl  +Gf         AzkAzk+i ,  lg. 

+    < 8  ^ZbCAzb+Azf)    +  e     8         ^ZbCAzb+Azf)      {0A*> 


AxwAyjAzk     n-i  IGel  -  Ge  Ax,2 

S=       XT UP         '      (  fA  )    Av    Av      ,     uee 


IGWI  +GW        Axj.i2 
"    (         16  AxwwAxiV    w 

jGni_1GIL    Ay,Ayi+i 

"    (        8  >  Aynn(Aynn  +  Ayn)  Unn 

,  IGsl  +  Gs  AyjAyu 

(  8  ;  Ayss(Ayss  +  Ays)  Uss 

(  IGfl  -  Gf  .        AzkAzk+1 
■    (         8         ;  Azff(Azff+  Azf)  Uff 
!Gbl  +  Gb  ,  AzkAzk,! 

(  8  }  Azbb(Azbb  +  Azb)  Ubb 

(PP-  Pw)AyjAzk  (6.19) 


AP=      AE  +    Aw+    AN+     As+     AF  +    Ab 

(    IGel  -  Ge  AXj2 

*         16        '  AxeAxi+i 

IGWI  +  GW        AxM2 
*         16  AxwwAxiV 

,   IGnl  -  Gn  Ay,Ayi+1 

<>  8  ;  Aynn(Aynn  +  Ayn) 

{  'Gsl  +  Gs  AyjAy;,, 

'    *  8  ;  Ayss(Ayss+ Ays) 

,  IGfl  -  Gf  AzkAzk+i 

1         8         }  Azff(Azff+  Azf) 

.   IGbl  +  Gb  v  AzkAzk.!  AxjAyjAzk  (f.  2m 

■     (  8  '  Azbb(Azbb  +  Azb)     +  At  ^ZU) 

6.1.3    v-Momentum  Equation  and  Linear  Equation  for  v 

The  coefficients  of  the  y-momentum  equations  are  the  following: 

Gn  pAxjAzk  IGnl  +  Gn      Ayi 

An  "    "    2     +    ^^yT  (         16         ;Ays 
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,  IGnl  -  Gn       Ay;  IG.I  -  Gs       Ay,!2 

+    (  8  }Ayj+i      +     (        16        >AysAyj 

Gs  pAxiAzk  ,  1G5I  +  Gs      Ay^ 

A*  =      T    +     PrA^T  +      (         8  }  AYj.2 

IG.I  -  Gs     AyM        fIGnl  +  Gn         Ayj2 

+    (        16        }Ays  +   <       16^      ;AysAyH 

GeAxj  pAysAzk  IGel  +  Ge  u AXlAxi+1 

Ae=    -^Ax7    +     Pl2Ax7    +     (  8  )(Axe(Axe  +  Axw) 

JGel-Ge    wAXiAxi+K  IGW1  -  Gw  AXjAXj.! 

+    < 8         )(  AxeAxee)    +     {  8         '  Axe(Axe  +  Axw) 


(6.21) 


(6.22) 


(6.23) 


GwAxj         PrAysAzk  1GWI  +  Gw        AxjAxm 

w  =        oav        +    n  9  A  v..       +      ^  8  M  Ax™Ax„™/ 


2Axv 


2Ax 


AxwAxv 


+    ( 


IGWI  -  Gw  ,        AxjAxj.i  IGel  +  Ge,        AxjAxi+i 


8 


) 


Axw(Axe+  Axw) 


8 


Axw(Axe  +  Axw) 


(6.24) 


Ap  ~    "    2Azf     +     n  2Azf 


+      ( 


IGfl  +  Gf 


AzkAzk+i 


IGfl  -  Gf    AzkAzk+i 
+    ^         8         ;  AzfAzff 


8  '  AzKAzb+  Azf) 

IGbl  -  Gb  v        AzkAzk.i 
^  8  ;  Azf(Azb  +  Azf) 


(6.25) 


AR  = 


2Azb 


+    t 


GbAzk 

2Azb 

IGbl  -  Gb  v       AzkAzk,i 
8         ;  Azb(Azb  +  Azf) 


IGbl  +  Gb  N  AzkAzk_i 


8 


AzbAzbb 


+   ( 


IGfl  +  Gf  N       AzkAzk+1 


8 


Azb(Azb  +  Azf) 


(6.26) 


S  = 


AxjAysAzk     n-i 


At 


vP 


IGnl  -  Gn 


16 
Gsl  +  Gs 


) 


Ayi 


AynAyJ+i 


VNN 


y 


j^L 


1  6         'AyssAyj  _2 
IGel  -  Ge  AxjAxi+i 

8  '    Axee(Axee+  Axe)       EE 

IGWI  +  Gw_  AxjAxj.j 

AXWw(/iXww  +  AxW/) 
AzkAzk+i 

)     K~.JK~.~l-    K~.\     VFF 


8 
IGfl  -  Gf 


8 

IGbl  +  Gb 
8 


Azff(Azff+  Azf) 


) 


AzkAzk.i 


Azbb(Azbb  +  Azb) 


vBb 


(PP  -  Ps  )A  XlA  zk  +  TsAyJAy^$AyHAxiAysAzk     (6-27) 


AP=      AE  +    Aw  +    AN  +     As+     AF  +    Ab 
IGnl  -  Gn  AVi2 

16         ;  AynAyJ+i 


( 
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2 


IG.I  +  Gs       Ayj.i"    , 
1        1 6  AynnAyj-2; 

IGe1  -  Ge  v  AxjAxj+i 

^  8  '    Axee(AXee  +  Axe) 

IGWI  +  Gw  .  AxjAxj.i 

^  8  ;  Axww(Axww  +  Axw) 

,  IGfi  -  Gf  AzkAzk+1 

(         8         ;  Azff(Azff+  Azf) 

.   IGhl  +  Gb  ,  AzkAzk.i  AxjAysAzk  6  2g) 

<  8  ;  Azbb(Azbb  +  Azb)  At 


6.1.4    z-Momentum  Equation  and  linear  equation  for  w 

The  coefficients  of  the  z-momentum  equations  are  the  following: 
a  Gf  pAxjAyj  IGfl+Gf      Agk. 

af  -  -  2    +   it  Azk     +    i      16      ;Azb 

IGfl  -  Gf      Azk  IGbl  -  Gb  v   Azk.!2  (629) 

+    (         8         }Az^  (         16         ^AzbAzk  V°-  yj 

Gb  pAxjAy;  IGbl  +  Gb  ,  Az^ 

A*  =      T   +     PrA^kT  +     (         8  }  Azk.2 

JGkJj^Gh  .Azk.i         ,IGnl  +  Gn         Azk2  (6  30) 

+    (         16        }Azb  +   (       16         'AzbAzn  K°-™} 


GeAxj          pAyjAzb  ,  'Gel  +  Ge  u Ax,Axi+1 

=    "  ~2Ax7    +     PlAx7~  C  8          M  Axe(Axe  +  Axw) 

.  IG,I  -  Ge       AxjAxi+i  IGwl  -  Gw     AxjAxj.i (6  3l) 

+    < 8 )(  AxeAxee  }  +  (          8          }  Axe(Axe  +  Axw)     ^U 


GwAxj         pAyjAzb  IGWI  +  Gw      ,  Ax.Ax^ 

Aw=      "2Ax7    +    ^lAxw     +     (  8  ^AxwAxwJ 

JGw<  -  Gw  AxjAxj,!  IGel  +  Ge ,        AxjAxi+i ,  ~2) 

+    *        8         }  Axw(Axe+  Axw)        * 8         ;Axw(Axe+Axw)         {0'J  } 

GnAyi  prAxjAzb  IGnl  +  Gn  ,  AyiAyi+i 

An=    "    2Ayn     +  ^^Ay^T  ("       8  >  Ayn(Ays+Ayn) 

lGnl  -  Gn  Ay,Ayi+i  ,      ,  'Gsl  -  Gs     AyjAy^ 

+    (  8  }AynA7n7  +      (         8         }  Ayn(Ayn  +  Ays)  ^j 

a  GsAy^         pAxjAzb  IG.I  +  Gs     AyiAyM 

As=       2Ays     +     ^   Ays       +      C        8  ;  AysAyss 

IG,I  -  Gs  AyjAyj.!  IGJ  +  Gn    AyiAyi+1 6  34 

+    i       8         >  Ays(Ays  +  Ayn)    +e    8         )  Ays(Ays  +  Ayn) 
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I 


AxjAyjAzb      n-i  IGfl  -  Gf         Azk2 

S"  At Wp  (        16        }AzfAzk+1  Wff 


.  IGbl  +  Gb  Azk.^ 

(~T6        ^'AzbbAzk+i  Wbb 
IGel  -  Ge  v  AxjAxj+i 

<  8  ^   Axee(AXee+  Axe)    WeE 

IGWI  +  GW  AxjAxj.! 

<>  8  ;  Axww(Axww  +  Axw)      w 

JGnL^Gn.     AyjAyi+1 

(  8         }  Aynn(Aynn+  Ayn)  Wnn 

<  8  '  Ayss(Ayss+  Ays)  Wss 

(PP-  PB)AxiAyj  (6.35) 


AP=      AE  +    Aw+    AN+     As+     AF  +    Ab 
(   IGfl  -  Gf  Azk2 

'         16        '  AzfAzk+i 
IGbL±Gb_  _AzM2_ 

'         16  AzffAzk-2 

IGel  -  Ge   s  AXjAXj+i 

*  8  '    Axee(Axee  +  Axe) 

IGWI  +  Gw  ,  AxjAxj.i 

*  8  '  Axww(Axww  +  Axw) 

,  IGnl  -  Gn  AyiAyi+1 

'    K         8  ;  Aynn(Aynn+  Ayn) 

IG,1  +  Gs  AyjAyj.!  AxjAyjAzb  (6  36) 

"     (  8  }  Ayss(Ayss+ Ays)     +  At  ^°J 


6.1.5    The  Pressure  Correction  Equation 

The  pressure  correction  equation  is  derived  from  the  discretized  continuity  equation. 
The  continuity  equation  is  the  following: 

du       dv        3w       n  ,,  o7n 

5-  +  3-    +3-    =0  (o.i/) 

dx       dy         dz 

Integrating  around  the  control  volume 

(  Ue  -  uw)A  yjA  zk    +       (  vn  -  vs)A  XiAzk     +       (  wf  -  wb)A  xtA  yj 
=     0  (6.38) 

A  pressure  correction  relationships  given  before  are  the  following: 
ue'  =due(PP' -  PE")      ; 
uw'  =duw(Pw' -  Pp)      ; 
vn   =dvn(PP' -  PN')       ; 
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vs'  =  dvs  OV  -  PP)  ; 
Wf'  =dwf(Pp'-  PF')  ; 
wb'  =  dwb(PB  -  PP')  (639) 

The  coupling  constants  are  calculated  using  the  SIMPLEX  algorithm.  Further  details  are 
given  in  Van  Doormal  and  Raithby  (1985).  also, 
ue  =  ue*  +  ue' 
uw  =  uw*  +  uw' 

Vn=    Vn*+    Vn 

vs  =  vs*  +  vs' 

Wf=    Wf*   +    Wf' 

wb  =  wb«  +  Wb1  <6-40> 

The  starred  quantities  represent  the  velocities  which  satisfy  the  momentum  equations  but 
not  the  continuity  equation.  The  primed  quantities  are  the  velocity  corrections  which  results 
in  velocities  that  satisfy  the  continuity  equation.  Hence,  using  equations  6.39  and  6.40  in 
6.38  we  obtain  a  pressure  correction  equation  of  the  following  type: 

ApPp  =  AEPE    +  AWPW    +  ANPN  +  ASPS'    +  AFPF    +  ABPB  +S 
Where, 

AE  =   dUeAyjAzk  ; 

Aw  =    duwAyjAzk  ; 

AN  =   dvnAx;Azk  ; 

As  =   dvsA  xjA  zk  ; 

AF  =   dwfA  xjA  yj  ; 

AB  =    dwbA  xjA  yj  ; 

AP=   AE    +    Aw    +    AN    +    As     +   AF    +    AB 

S  =     -ue*AyjAzk  +   uw*AyjAzk  -   vn*A  x^  zk  +  vs*Ax,Azk 

-   wf*A  xjA  yj  +   wb*A  x;A  yj  (6.41) 

After  solving  the  equations  for  pressure  correction,  the  velocities  and  pressures  are 
corrected. 

6.2       Numerical  Boundary  Conditions 

All  the  coefficients  and  the  relations  shown  so  far  are  valid  for  interior  points.  Near 
the  boundary  some  special  treatment  is  necessary.  For  example,  the  central  differencing 
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used  to  estimate  the  diffusive  fluxes  cannot  be  used  for  the  wall  since  there  is  no  grid  point 
beyond  the  wall.  There  are  basically  two  approaches  open  to  us: 

1 .  Use  different  interpolation  and  differencing  scheme  as  we  approach  the  wall.  For 
instance,  a  one  sided  three  point  differencing  scheme  could  be  used  at  the  wall. 

2.  Use  the  same  interpolation  and  differencing  scheme  but  create  additional  fictitious 
grid  points  or  control  volumes  beyond  the  boundary  and  assign  suitable  values  to  the  grid 
points  for  the  different  variables.  Both  approaches  are  equivalent. 

The  second  approach  was  used  for  this  problem.  The  advantage  being  that  it  is  a  lot 
simpler  to  implement  since  the  same  scheme  is  used  with  modifications.  To  that  end, 
additional  control  volumes  were  used  outside  the  calculation  domain  to  impose  the 
boundary  condition.  The  one  guideline  we  follow  is  that  the  boundary  conditions  imposed 
should  be  second  order  or  higher. 

6.2.1      Velocity  Boundary  Conditions 

The  no  slip  boundary  conditions  are  to  be  imposed  at  all  walls.  For  the  x- 
momentum  equation,  the  u-velocity  grid  is  staggered  in  the  x-direction  the  grid  point 
coincides  with  the  wall  for  the  east  and  west  walls  perpendicular  to  the  x-direction  as 
shown  in  figure.  In  the  other  directions  the  control  volume  surfaces  coincide  with  the  walls 
as  shown  in  figure.  For  the  staggered  direction  the  computational  grid  extends  up  to  the 
point  i  as  shown.  Point  i+1  is  the  wall,  i+2  is  the  pseudo-point  beyond  the  wall.  Up  to  i-1 
no  special  treatment  is  necessary.  For  the  equation  of  point  i  a  special  value  needs  to  be 
assigned  to  i+2.  Noting  that  the  wall  is  a  plane  of  anti-symmetry  for  a  scheme  of  second 
order  accuracy  we  impose, 

Ui+2  =  "Ui 

m+i  is  of  course  set  equal  to  zero.  Consequently,  substituting  the  values  in  equation  A-20 
for  the  point  i , 

AE  =  0,  the  source  term  S  is  affected  also  because  uEE  =  -uP .  Similarly,  AP  is  suitably 
modified.  The  west  wall  is  similarly  dealt  with. 

For  the  other  walls  only  the  control  volume  next  to  the  wall  needs  to  be  dealt  with 
hence  ui+i  =  -uj  or, 
Up  =  -  uE 

therefore,  AP  =    AP  +  AE  and  AE  =0 

The  v-velocity  and  the  w- velocity  boundary  conditions  are  similarly  dealt  with.  What  must 
be  kept  in  mind  is  that  the  v-velocity  grid  is  staggered  in  the  y  direction  and  that  the  w- 
velocity  grid  is  staggered  in  the  z  direction. 
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In  cases  where  symmetry  was  exploited  uj+i  =  uj  was  used  instead. 

6.2.2  Temperature  Boundary  Conditions 

There  are  two  types  of  boundary  conditions  to  be  considered;  adiabatic  and 
isothermal.  The  temperature  grid  is  not  staggered.  The  adiabatic  boundary  condition  is 
identical  to  the  symmetry  condition  hence  for  the  grid  point  next  to  the  east  wall 
TP  =TE 

therefore,  AP  =   AP  -  AE  and  AE  =  0. 

For  the  isothermal  boundary  condition  the  value  of  the  temperature  grid  point  outside  the 
domain  was  assigned  by  a  linear  extrapolation.  Hence,  we  have  the  following; 

TE  =    2Twaii  -  lp 
The  symmetry  conditions  are  identical  to  the  adiabatic  conditions. 

6.2.3  Pressure  Correction  Boundary  Conditions 

The  velocity  corrections  at  the  wall  must  all  be  zero.  This  is  because  as  a  result  of 
the  no-slip  conditions  the  velocities  at  the  wall  are  identically  zero.  Hence, 

Pp=Pe' 

Therefore,  in  the  governing  equations; 

AP  =    AP  -  AE  and  AE  =  0 
Similar  modifications  will  have  to  be  made  at  the  other  walls. 
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APPENDIX   C 


6. 1       General  Remarks 

The  listing  of  the  code  is  provided  in  this  appendix.  The  code  is  written  in  a  manner 
that  facilitates  comprehension.  The  indentation  of  all  DO  loops  and  space  between  operands 
all  add  to  the  clarity.  There  are  plenty  of  comments  inserted  in  the  code.  In  fact,  it  is  quite 
possible  to  understand  most  of  the  code,  without  actually  going  through  the  documentation. 

Note  that  all  three-dimensional  variables  are  subscripted  (NIP1,NJP1,NKP1).  A 
higher  number  can  also  be  dimensioned,  but  that  is  a  waste  of  computer  resources.  It  is 
assumed  that  the  user  has  access  to  some  full-screen  editors,  therefore  it  is  quite 
straightforward  to  change  indices  for  the  arrays.  That  may  be  required  when  the  grid  sizes 
are  changed  or  if  a  different  physical  problem  is  being  attempted.  All  one-dimensional 
variables  in  the  subroutine  SIP  must  be  dimensioned  NNMAX  which  is  a  product  of  NIP1, 
NJP1  and  NKP1.  In  the  specific  example  of  the  listing,  NIP1  =  12,  NJP1  =  16,  NKP1  = 
16  and  NNMAX  =  3072.  This  problem  requires  about  1.2  Mbytes  of  memory  in  a 
RS/6000  workstation.  The  listing  follows. 
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PROGRAM  AMPHIB 
c 

C 

CALL  OPENF 
c 

call  grid 
c 

call  prop 
c 

call  initio 
c 

call  ploop 
c 

c*****^*************************************************************** 
BLOCK  DATA 

implicit  real*8  (a-h,o-z) 
common/blayer /xbr , ybr , zbr 
common/unf rm/ iunf rm 

common/parm/ra , pr , sorsum , dtime , xper , roll , dt_inv , xt lme 
common/tol/small , eps, sormax 
common/dims/h,wth,bth,hchip,wchip,bchip,bsub 

common/count/nt,nmax, imax, itmax, krun, nprint 

common/chip/ ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) ,nchp, 
&  isub, ichip,jchip,kchip 

common/ limits/ni , nipl , niml , nj , n j pi, njml , nk, nkpl, nkml , 
&  nip2,njp2,nkp2, iter,nnmax 

common/diff/alc,als,thot,tcool,tavg 

common/scheme/quick, upwind 

common/array/njchip, nkchip, ichoice 

common/spaced/ychip(3) ,zchip(3) 

COMMON/RHOCP/RHS , RHC 

COMMON/ POWER/QQQ,QCOND 

c 

DATA  NIP2,NIP1,NI,NIM1/13, 12,11,10/ 

DATA  NJP2,NJP1,NJ,NJM1/17, 16, 15,14/ 
DATA  NKP2 , NKP1 , NK , NKM1/ 17 ,16,15,14/ 

DATA  NNMAX , NMAX , IMAX , ITMAX , KRUN , NPRINT/ 3072,4000,10,5, 
&  1,100/ 

DATA  SMALL , EPS , SORMAX/ 1 . 0e2 0 , 1 . OE-8 , 2 . / 

DATA  ISUB,ICHIP,JCHIP,KCHIP,NCHP,ICHOICE/2,2,2,2,l,l/ 
DATA  XBR,YBR,ZBR/1. 50, 1.00, 1.50/ 
DATA  H,WTH,BTH/140. ,140.0, 100./ 
DATA  HCHIP,WCHIP,BCHIP,BSUB/24. ,8. ,6. ,19.5/ 
DATA  ALC,ALS,THOT,TCOOL,TAVG/3338. , 2 . 0 , 0 . 0 , 0 . 0 , 19 . 0/ 
DATA  RHC,RHS/1. 20, 0.262/ 
DATA  QQQ , DTIME , pr , ra , XPER, ROLL/ 2 . 0e-03 , 2 . 0E-7 , 28 . , 1 . 0e6 , 0 . 0 , 0 . 0/ 

DATA  IUNFRM/1/ 

DATA  QUICK/ 1./ 

DATA  NJCHIP,NKCHIP/3,3/ 

DATA  YCHIP,ZCHIP/3  8.0,7  6.0,114.0,50.8,101.6,152.4/ 

END 
c  ********************************************************************* 

SUBROUTINE  CALT 

C 

C      THIS  SUBROUTINE  ASSEMBLES  THE  DISCRETE  ANALOG  OF  THE  ENERGY 

C      EQUATION,  APPLIES  THE  BOUNDARY  CONDITION  AND  THEN  SOLVES  FOR 

C      THE  TEMPERATURE. 

C 
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c 
c 

c 


implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ra,pr,sorsum,dtime,xper,roll,dt_inv,xtime 

common/ limits/ni , nipl , niml , nj , njpl , njml , nk, nkpl , nkml , 
&  nip2 , njp2 , nkp2 , iter , nnmax 

common/fv_init/tod(12,16,16) ,uod ( 12 , 16 , 16) , vod ( 12 , 16, 16) , 
&  wod(12,16,16) ,pod(12,16,16) 

common/fv_cur/t(12,16,16) ,u(12,16,16),v(12,16,16), 
&  w(12,16,16) ,p(12,16,16) 

common/fv_int/tpd(12,16,16) , upd ( 12 , 16 , 16) , vpd ( 12 , 16 , 16) , 
&  wpd(12, 16, 16) ,ppd(12, 16, 16) 

common/coeff/ap(12, 16, 16) ,aw(12,16,16) ,ae(12,16,16) , 
&  as (12, 16, 16) ,an(12,16,16) ,ab(12,16, 16) , 

&  af  (12,16,16)  ,su(12, 16, 16) 

COMMON /COEF2/AWW( 12, 16, 16) ,AEE(12,16, 16) , ASS (12, 16, 16) , 
&  ANN(12,16,16) ,ABB(12,16,16) , AFF ( 12 , 16 , 16) 

common/chip/ ibgn,iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) , 
&  isub, ichip, jchip,kchip 

common/condu/alpha(12, 16, 16) ,rho(12 

COMMON / POWER / QQQ , QCOND 

COMMON/SCHEME/QUICK, UPWIND 

common/dif f /ale, als, thot , tcool, tavg 

common/array/njehip, nkchip, ichoice 

CALCULATE  COEFFICIENTS 

do  80  k=2,nk 
do  80  j=2,nj 

do  80  1=2, ni 

VOLDT  =  DXX(I)  *  DYY(J)  *  DZZ(K)  /  DTIME 
&  *  RHO(I,J,K) 


nchp, 
16, 16) , visco(12, 16, 16) 


c 
c 

CALCULATE  INTERFACIAL  THERMAL  DIFFUSIVITIES 

c 

USING  THE 

HARMONIC  MEAN  FORMULATION 

c 

CONDW 

- 

ALPHA (I, J, K)  *  ALPHA(I-1,J,K) 

& 

* 

(DXX(I-l)  +  DXX(I))  /  (DXX(I-l) 

& 

* 

ALPHA(I,J,K)  +  DXX(I)  *  ALPHA ( 1-1 , J , K) ) 

CONDE 

= 

ALPHA(I,J,K)  *  ALPHA(I+1,J,K) 

& 

* 

(DXX(I+1)  +DXX(I))  /  (DXX(I+1) 

& 

* 

ALPHA(I,J,K)  +  DXX(I)  *  ALPHA ( 1+1 , J , K) ) 

CONDS 

= 

ALPHA(I,J,K)  *  ALPHA(I, J-1,K) 

& 

* 

(DYY(J-l)  +  DYY(J))  /  (DYY(J-l) 

& 

* 

ALPHA(I,J,K)  +  DYY(J)  *  ALPHA ( I , J-l , K) ) 

CONDN 

= 

ALPHA(I,J,K)  *  ALPHA(I,J+1,K) 

& 

* 

(DYY(J+1)  +  DYY(J))  /  (DYY(J+1) 

& 

* 

ALPHA(I,J,K)  +  DYY(J)  *  ALPHA ( I , J+l , K) ) 

CONDB 

= 

ALPHA(I,J,K)  *  ALPHA(I, J,K-1) 

& 

* 

(DZZ(K-l)  +  DZZ(K))  /  (DZZ(K-l) 

& 

* 

ALPHA(I,J,K)  +  DZZ(K)  *  ALPHA ( I , J , K-l) ) 

CONDF 

= 

ALPHA(I,J,K)  *  ALPHA(I,J,K+1) 

& 

* 

(DZZ(K+1)  +  DZZ(K))  /  (DZZ(K+1) 

& 

* 

ALPHA(I,J,K)  +  DZZ(K)  *  ALPHA ( I , J , K+l) ) 

c 

CONDN1 

=  DZZ(K)  *  DXX(I)  /  DYYS(J+1)  *  CONDN 

CONDS1 

=  DZZ(K)  *  DXX(I)  /  DYYS(J)  *  CONDS 

CONDF1 

=  DXX(I)  *  DYY(J)  /  DZZS(K+1)  *  CONDF 

CONDB1 

=  DXX(I)  *  DYY(J)  /  DZZS(K)  *  CONDB 
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C0NDE1  =  DYY(J)  *  DZZ(K)  /  DXXS(I+1)  *  CONUfc. 
CONDW1  =  DYY(J)  *  DZZ(K)  /  DXXS(I)  *  CONDW 

GP  =  RHO(I,J,K) 

GW  =  RHO(I-l,J,K) 

GE  =  RHO(I+l,J,K) 

GS  =  RHO(I,J-l,K) 

GN  =  RHO(I,J+l,K) 

GB  =  RHO(I,J,K-l) 

GF  =  RHO(I,J,K+l) 

RW  =  (GP  *  DXX(I-l)  +  GW  *  DXX(I)) 

&  /  (DXX(I-l)  +  DXX(I)) 

RE  =  (GP  *  DXX(I+1)  +  GE  *  DXX(I)) 

&  /  (DXX(I+1)  +  DXX(I) ) 

RS  =  (GP  *  DYY(J-l)  +  GS  *  DYY(J)) 

&  /  (DYY(J-l)  +  DYY(J)) 

RN  =  (GP  *  DYY(J+1)  +  GN  *  DYY(J)) 

&  /  (DYY(J+1)  +  DYY(J) ) 

RB  =  (GP  *  DZZ(K-l)  +  GB  *  DZZ(K)) 

&  /  (DZZ(K-l)  +  DZZ(K) ) 

RF  -  (GP  *  DZZ(K+1)  +  GF  *  DZZ(K)) 

&  /  (DZZ(K+1)  +  DZZ(K) ) 

CE  =  U(I+1,J,K)  *  DYY(J)  *  DZZ(K)  *  RE 

CW  =  U(I,J,K)  *  DYY(J)  *  DZZ(K)  *  RW 

CN  =  V(I,J+1,K)  *  DZZ(K)  *  DXX(I)  *  RN 

CS  =  V(I,J,K)  *  DZZ(K)  *  DXX(I)  *  RS 

CF  =  W(I,J,K+1)  *  DXX(I)  *  DYY(J)  *  RF 

CB  =  W(I,J,K)  *  DXX(I)  *  DYY(J)  *  RB 

AE(I,J,K)  =  CONDE1  +  QUICK 

&  *  (-  0.5  *  CE  *  DXX(I)  /  DXXS(I+1) 

&  +  0.125  *  (abs(ce)  +  ce)  *  dxx(i)  *  dxx(i+l) 

&  /  (dxxs(i+l)  *  (dxxs(i)  +dxxs(i+l))) 

&  +  0.125  *  (abs(ce)  -  ce)  *  dxx(i)  *  dxx(i+l) 

&  /  (dxxs(i+l)  *  dxxs(i+2)) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (DXXS(I+1)  *  (DXXS(I)  +  DXXS(I+1)))) 

&  +  UPWIND  *  (  0.5  *  (ABS(CE)  -  CE) ) 

AW(I,J,K)  =  CONDW1  +  QUICK 

&  *  (0.5  *  CW  *  DXX(I)  /  DXXS(I) 

&  +  0.125  *  (abs(cw)  +  cw)  *  dxx(i)  *  dxx(i-l) 

&  /  (dxxs(i-l)  *  dxxs(i) ) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (dxxs(i)  *  (dxxs(i)  +  dxxs(i+l))) 

&  +  0.125  *  (abs(ce)  +  ce)  *  dxx(i+l)  *  dxx(i) 

&  /  (DXXS(I)  *  (DXXS(I)  +  DXXS(I+1)))) 

&  +  UPWIND  *  (  0.5  *  (ABS(CW)  +  CW) ) 

AN(I,J,K)  =  CONDN1  +  QUICK 

&  *  (-  0.5  *  CN  *  DYY(J)  /  DYYS(J+1) 

&  +  0.125  *  (abs(cn)  +  en)  *  dyy(j)  *  dyy(j+l) 

&  /  (dyys(j+l)  *  (dyys(j)  +  dyys(j+l))) 

&  +  0.125  *  (abs(cn)  -  en)  *  dyy ( j )  *  dyy(j+l) 

&  /  (dyys(j+l)  *  dyys(j+2)) 

&  +  0.125  *  (abs(cs)  -  cs)  *  dyy(j-l)  *  dyy ( ] ) 

&  /  (DYYS(J+1)  *  (DYYS(J)  +  DYYS(J+1)))) 

&  +  UPWIND  *  (  0.5  *  (ABS(CN)  -  CN) ) 

AS(I,J,K)  =  CONDS1  +  QUICK 
&  *  (0.5  *  CS  *  DYY(J)  /  DYYS(J) 
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+  0.125  *  (abs(cs) 
/  (dyys(j-l)  *  dyy 
+  0.125  *  (abs(cs) 
/  (dyys(j)  *  (dyys 
+  0.125  *  (abs(cn) 
/  (DYYS(J)  *  (DYYS 
+  UPWIND  *  (  0.5  * 
AF(I,J,K)  =  CONDF1  +  QUICK 

*  (-  0.5  *  CF  *  DZ 

0.125  *  (abs(cf) 

(dzzs(k+l)  *  (dz 

0.125  *  (abs(cf) 

(dzzs(k+l)  *  dzz 

0.125  *  (abs(cb) 

/  (DZZS(K+1)  *  (DZ 

+  UPWIND  *  (  0.5  * 

AB(I,J,K)   >  CONDB1  +  QUICK 

(0.5  *  CB  *  DZZ( 

0.125  *  (abs(cb) 

(dzzs(k-l)  *  dzz 

0.125  *  (abs(cb) 

/  (dzzs(k)  *  (dzzs 

+  0.125  *  (abs(cf) 

/  (DZZS(K)  *  (DZZS 

+  UPWIND  *  (  0.5  * 


+ 

/ 

+ 

/ 

+ 


+ 

/ 

+ 


+  cs)  *  dyy(j)  *  dyyCj-l) 

s(j)) 

-  cs)  *  dyy(j-l)  *  dyy(j) 
(j)  +  dyys(j+l) ) ) 

+  en)  *  dyy(j+l)  *  dyy ( j ) 
(J)  +  DYYS(J+1) ) ) ) 
(ABS(CS)  +  CS) ) 

Z(K)  /  DZZS(K+1) 

+  cf)  *  dzz(k)  *  dzz(k+l) 
zs(k)  +  dzzs(k+l) ) ) 

-  cf)  *  dzz(k)  *  dzz(k+l) 
s(k+2) ) 

-  cb)  *  dzz(k-l)  *  dzz(k) 
ZS(K)  +  DZZS(K+1) ) ) ) 

(ABS(CF)  -  CF) ) 

K)  /  DZZS(K) 

+  cb)  *  dzz(k)  *  dzz(k-l) 
s(k)) 

-  cb)  *  dzz(k-l)  *  dzz(k) 
(k)  +  dzzs(k+i) ) ) 

+  cf)  *  dzz(k+l)  *  dzz(k) 
(K)  +  DZZS(K+1) ) ) ) 
(ABS(CB)  +  CB) ) 


AEE(I, J,K) 


C 
C 

c 


/ 

AWW(I,J,K)  = 
* 

/ 
ANN(I,J,K)  = 


ASS(I,J,K)  = 
* 

/ 
AFF(I,J,K)  = 

/ 

ABB(I,J,K)  = 

* 

/ 


& 
& 

& 
& 

& 
& 

& 
& 

& 
& 
i 
& 
& 


-  0.125  *  QUICK 
(ABS(CE)  -  CE)  *  DXX(I) 
(DXXS(I+2)  *  (DXXS(I+1) 

-  0. 125  *  QUICK 
(ABS(CW)  +  CW)  *  DXX(I) 
(DXXS(I-l)  *  (DXXS(I-l) 

-  0.125  *  QUICK 
(ABS(CN)  -  CN)  *  DYY(J) 
(DYYS(J+2)  *  (DYYS(J+1) 

-  0.125  *  QUICK 
(ABS(CS)  +  CS)  *  DYY(J) 
(DYYS(J-l)  *  (DYYS(J-l) 

-  0.125  *  QUICK 
(ABS(CF)'  -  CF)  *  DZZ(K) 
(DZZS(K+2)  *  (DZZS(K+1) 

-  0. 125  *  QUICK 
(ABS(CB)  +  CB)  *  DZZ(K) 
(DZZS(K-l)  *  (DZZS(K-l) 


ap(i,j,k)  =  ae(i,j,k)  +  aw(i,j,k)  + 

&  +  AF(I,J,K)  +  AB(I,J,K)  + 

&  +  AEE(I,J,K)  +  ASS(I,J,K) 

&  +  ABB(I,J,K)  +  AFF(I,J,K) 

SU(I,J,K)  =  VOLDT  *  TOD (I, J, K) 

80   continue 

CALCULATION  OF  THE  SOURCE  TERM 


DO  81  1=4, NI 


DO 

81 

J=2 

,NJ 

DO 

81 

K=2 

NK 

SU( 

I, J 

K) 

CONTINUE 

*  DXX(I+1) 

+  DXXS(I+2) ) ) 

*  DXX(I-l) 
+  DXXS(I) ) ) 

*  DYY(J+1) 

+  DYYS(J+2) ) ) 

*  DYY(J-l) 
+  DYYS(J) ) ) 

*  DZZ(K+1) 

+  DZZS(K+2) ) ) 

*  DZZ(K-l) 
+  DZZS(K) ) ) 


an(i,j,k)  +  as(i, j ,k) 
AWW(I,J,K) 

+  ANN(I,J,K) 
+  VOLDT 


=  SU(I,J,K) 
+  AWW(I,J,K) 


*  TPD(I-2,J,K) 


63 


DO  82  I=2,NI-2 
DO  8  2  J=2,NJ 

DO  82  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AEE(I,J,K)  *  TPD(I+2,J,K) 

82  CONTINUE 
C 

DO  83  1=2, NI 

DO  8  3  J=4,NJ 

DO  83  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ASS(I,J,K)  *  TPD(I,J-2,K) 

83  CONTINUE 
C 

DO  84  1=2, NI 

DO  84  J=2,NJ-2 
DO  84  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ANN(I,J,K)  *  TPD(I,J+2,K) 

84  CONTINUE 
C 

DO  85  1=2, NI 

DO  8  5  J=2,NJ 

DO  85  K=4,NK 

SU(I,J,K)  -  SU(I,J,K) 
&  +  ABB(I,J,K)  *  TPD(I,J,K-2) 

85  CONTINUE 
C 

DO  86  1=2, NI 

DO  8  6  J=2,NJ 

DO  86  K=2,NK-2 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AFF(I,J,K)  *  TPD(I,J,K+2) 

86  CONTINUE 
C 

if  (nchp  .ne.  0)  then 
C 

C        CHIP  HEAT  GENERATION 
C 

do  90  m=l,njchip 
do  90  n=l,nkchip 

do  90  i=ibgn,iend 

do  90  j=jbgn(m) , jend(m) 

do  90  k=kbgn(n) ,kend(n) 

su(i,j,k)  =  ale  *  dxx(i)  *  dyy ( j ) 
&  *  dzz (k)  +  su(i, j ,k) 

90       continue 
endif 
C 
C         BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM, ISOTHERMAL) 

c 

do  101  k=2,nk 
do  101  i=2,ni 

AP(I,2,K)  =  AP(I,2,K)  +  AS(I,2,K)  -  ASS(I,2,K) 
su(i,2,k)  =  su(i,2,k)  +  2.  *  as(i,2,k)  *  thot 
as(i, 2 ,k)  =  0. 

SU(I,3,K)  =  SU(I,3,K)  +  ASS(I,3,K) 
&  *  (2.  *  THOT  -  TPD(I,2,K)) 

SU(I,NJM1,K)  =  SU(I,NJM1,K)  +  ANN ( I , NJM1 , K) 
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&  *  (2.  *  TCOOL  -  TPD(I,NJ,K)) 

AP(I,NJ,K)  =  AP(I,NJ,K)  +  AN(I,NJ,K)  -  ANN(I,NJ/K) 
su(i,nj,k)  =  su(i,nj,k)  +  2.  *  an(i,nj,k)  *  tcool 
an(i, nj ,k)  =  0. 

101  continue 

c 

c         ADIABATIC  WALLS  (LEFT  AND  RIGHT) 

c 

do  102  k=2,nk 
do  102  j=2,nj 

AP(2,J,K)  =  AP(2,J,K)  -  AW(2,J,K)  -  AWW(2,J,K) 
aw(2,j,k)  =  0. 

SU(3,J,K)  =  SU(3,J,K)  +  AWW(3,J,K)  *  TPD(2,J,K) 
AP(NI,J,K)  =  AP(NI,J,K)  -  AE(NI,J,K)  -  AEE(NI,J,K) 
SU(NIM1,J,K)  =  SU(NIM1,J,K) 
&  +  AEE(NIM1,J,K)  *  TPD(NI,J,K) 

ae(ni, j,k)  =  0. 

102  continue 

c 

c         ADIABATIC  WALLS  (BACK  AND  FRONT) 

c 

do  103  j=2,nj 

do  103  i=2,ni 

AP(I,Jf2)  =  AP(I,J,2)  -  AB(I,J,2)  -  ABB (I, J, 2) 

ab(i,j,2)  =  0. 

SU(I,J,3)  =  SU(I,J,3)  +  ABB(I,J,3)  *  TPD(I,J,2) 
AP(I,J,NK)  =  AP(I,J,NK)  -  AF(I,J,NK)  -  AFF(I,J,NK) 
SU(I,J,NKM1)  =  SU(I,J,NKM1) 
&  +  AFF(I,J,NKM1)  *  TPD(I,J,NK) 

af(i,j,nk)  =  0. 

103  continue 
c 

C        SOLVE  DIFFERENCE  EQUATION 

C 


call  sip  (2,2, 2,nifnj,nk,t) 


c 


C      PRESCRIBE  WALL  TEMPERATURES  CONSISTENT  WITH  THE 
C  BOUNDARY  CONDITIONS 

C  (RIGHT  AND  LEFT  WALL) 

C 

do  310  k=2,nk 
do  310  j=2,nj 

t(nipl, j,k)  =  t(ni, j ,k) 
t(l,j,k)  =  t(2,j,k) 
310  continue 
c 

c  FRONT  AND  BACK  WALL 

c 

do  320  j=2,nj 

do  320  i=2,ni 

t(i, j,nkpl)  =  t(i, j ,nk) 
t(i,j,l)  =  t(i,j,2) 
320  continue 
c 

c  TOP  AND  BOTTOM  WALL 

c 

do  330  k=2,nk 

do  330  i=2,ni 

t(i,l,k)  =  2.  *  thot  -  t(i,2,k) 
t(i,njpl,k)  =  2.  *  tcool  -  t(i,nj,k) 


65 


330  continue 

c      THE  KINEMATIC  VISCOSITY  IS  CALCULATED,  WHICH  IS  A  STRONG 
c      FUNCTION  OF  TEMPERATURE 

c 

do  350  i=2,ni 
do  350  j=2,nj 

do  350  k=2,nk 

ttmp  =  tavg  +  qcond  *  t(i,j,k) 
if  (ichoice  .eq.  0)  then 

gnu  =  pr 
endif 
if  (ichoice  .eq.  1)  then 

gnu  =  (1.4074  -  2.96e-2  *  ttmp 
&  +  3.8018e-4  *  ttmp**2 

&  -  2.731e-6  *  ttmp**3 

5  +  8.168e-9  *  ttmp**4)  *  29.088 

endif 

if  (ichoice  .eq.  2)  then 

gnu  =  (251.62  -  13.723  *  tavg 

6  +  3.056e-l  *  tavg**2 

&  -  3.1704e-3  *  tavg**3 

+  1.2668e-5  *  tavg**4)  *  28.666 


& 


350  continue 


endif 
visco(i,j,k)  =  gnu 


SUBROUTINE  CALU 

C 

C      THIS  SUBROUTINE  ASSEMBLES  THE  U-MOMENTUM  EQUATION  AND 

C      BOUNDARY  CONDITION  AND  SOLVES  FOR  THE  U-VELOCITY  (X-DIRECTION) 

C 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ra , pr , sorsuro , dt  ime , xper , roll , dt_inv , xtime 

common/ limits/ni , nipl , niml , nj , njpl, njml , nk, nkpl, nkml, 
&  nip2,njp2,nkp2 , iter,nnmax 

common/fv_init/tod(12,16,16) , uod ( 12 , 16 , 16) , vod ( 12 , 16 , 16) , 
&  wod(12,16,16) ,pod(12,16,16) 

common/fv_cur/t(12,16,16) ,u(12,16,16),v(12,16,16), . 
&  w(12,16f 16) ,p(12, 16,16) 

common/fv_int/tpd(12,16,16) , upd ( 12 , 16 , 16) , vpd ( 12 , 16 , 16) , 
&  wpd(12,16,16) ,ppd(12,16, 16) 

common/mscn/smp(12, 16, 16) ,resorm(93) , 
&  du(12,16, 16) ,dv(12, 16,16) , 

&  dw(12, 16, 16) ,pp(12, 16, 16) 

common/coeff/ap(12, 16, 16)  ,aw(12,16,16)  ,ae(12,16,16)  , 
&  as ( 12, 16, 16 ), an (12, 16,16) ,ab( 12, 16,16) , 

&  af (12, 16, 16) ,su(12, 16, 16) 

COMMON/COEF2/AWW(12, 16, 16)  ,AEE(12, 16, 16)  , ASS (12, 16, 16)  , 
&  ANN(12,16#16) ,ABB(12,16,16) , AFF ( 12 , 16 , 16) 

common/condu/alpha(12,16,16) , rho( 12 , 16, 16) , visco ( 12 , 16 , 16) 

common/chip/ ibgn, iend, jbgn(3) , jend(3) , kbgn ( 3 ) ,kend(3) ,nchp, 
&  isub, ichip, jchip,kchip 

common/scheme/quick, upwind 
common/dif f /ale, a Is, thot , tcool , tavg 
common/array/njehip, nkchip, ichoice 
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V1SW  = 

vise  = 

viss  = 

& 

* 

& 

* 

visn  = 

& 

* 

& 

* 

visb  = 

& 

* 

& 

* 

visf  = 

& 

* 

& 

* 

common/tol/small,eps,sormax 
common/ if irst/nust , nvst , nwst , npst 

c 

c         CALCULATE  COEFFICIENTS 

c 

do  100  k=2,nk 
do  100  j=2,nj 

do  100  i=nust,ni 

voldt  =  dxxs(i)  *  dyy(j)  *  dzz(k)  /  dtime 

C 

C        CALCULATE  INTERFACIAL  KINEMATIC  VISCOSITIES 

C  USING  THE  HARMONIC  MEAN  FORMULATION 

C 

visco(i-l, j ,k) 

visco( i , j ,k) 

visco(i,j,k)  *  visco (i,j-l,k) 
(dyy(j-l)  +dyy(j))  /  (dyy(j-l) 
visco(i,j,k)  +  dyy(j)  *  visco( l , j-1 , k) ) 
visco(i, j,k)  *  visco(i, j+l,k) 
(dyy(j+l)  +dyy(j))  /  (dyy(j+l) 
visco(i,j,k)  +  dyy(j)  *  visco ( l , j+1 , k) ) 
visco(i,j,k)  *  visco(i, j ,k-l) 
(dzz(k-l)  +dzz(k))  /  (dzz(k-l) 
visco(i,j,k)  +  dzz(k)  *  visco ( i , j , k-1) ) 
visco(i, j ,k)  *  visco(i, j ,k+i) 
(dzz(k+l)  +dzz(k))  /  (dzz(k+l) 
visco(i,j,k)  +  dzz(k)  *  visco ( i , j , k+i) ) 

dxys  =  dxxs(i)  *  dyy(j) 
dyzs  =  dyy(j)  *  dzz(k) 
dzxs  =  dzz(k)  *  dxxs(i) 
zxoyn  =  dzxs  /  dyys(j+l) 
zxoys  =  dzxs  /  dyys(j) 
xyozf  =  dxys  /  dzzs(k+l) 
xyozb  =  dxys  /  dzzs(k) 
yzoxe  =  dyzs  /  dxx(i) 
yzoxw  =  dyzs  /  dxx(i-l) 

GN  =  V(I,J+1,K) 
GNW  =  V(I-1,J+1,K) 
GS  =  V(I,J,K) 
GSW  =  V(I-1,J,K) 

GE  =  U(I+1,J,K) 
GP  =  U(I,J,K) 
GW  =  U(I-1,J,K) 

GF  =  W(I,J,K+1) 
GFW  =  W(I-1,J,K+1) 
GB  =  W(I,J,K) 
GBW  =  W(I-1,J,K) 

en  =  (gn  *  dxx(i-l)  +  gnw  *  dxx(i)) 
&  /  (dxx(i-l)  +  dxx(i))  *  dzxs 

cs  =  (gs  *  dxx(i-l)  +  gsw  *  dxx(i)) 
&  /  (dxx(i-l)  +  dxx(i))  *  dzxs 

ce  =  0.5  *  (ge  +  gp)  *  dyzs 

cw  =  0.5  *  (gp  +  gw)  *  dyzs 

cf  =  (gf  *  dxx(i-l)  +  gfw  *  dxx(i)) 
&  /  (dxx(i-l)  +  dxx(i))  *  dxys 
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cb  =  (gb  *  dxx(i-l)  +  gbw  *  dxx(i)) 
&  /  (dxx(i-l)  +  dxx(i))  *  dxys 

visel  =  yzoxe  *  vise 

viswl  =  yzoxw  *  visw 

visnl  =  zxoyn  *  visn 

vissl  =  zxoys  *  viss 

visfl  =  xyozf  *  visf 

visbl  =  xyozb  *  visb 

AE(I,J,K)  =  VISE1  +  (-  0.5  *  CE  +  0.0625  *  (ABS(CE)  +  CE) 

&  *  dxx(i)  /  dxxs(i)  +  0.125  *  (abs(ce)  -  ce) 

&  *  dxx(i)  /  dxx(i+l)  +  0.0625  *  (abs(cw)  -  cw) 

&  *  DXX(I-1)**2  /  (DXX(I)  *  DXXS(I)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CE)  -  CE) ) 

AW(I,J,K)  =  VISW1  +  (0.5  *  CW  +  0.125  *  (ABS(CW)  +  CW) 

&  *  dxx(i-l)  /  dxx(i-2)  +  0.0625  *  (abs(cw)  -  cw) 

&  *  dxx(i-l)  /  dxxs(i)  +  0.0625  *  (abs(ce)  +  ce) 

&  *  DXX(I)**2  /  (DXX(I-l)  *  DXXS(I)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CW)  +  CW) ) 

AN(I,J,K)  =  VISN1  +  (-  0.5  *  CN  *  DYY(J)  /  DYYS(J+1) 

&  +  0.125  *  (abs(cn)  +  en)  *  dyy ( j )  *  dyy(j+l) 

&  /  (dyys(j+l)  *  (dyys(j)  +dyys(j+l))) 

&  +  0.125  *  (abs(cn)  -  en)  *  dyy ( j )  *  dyy(j+l) 

&  /  (dyys(j+l)  *  dyys(j+2) ) 

&  +  0.125  *  (abs(cs)  -  cs)  *  dyy(j-l)  *  dyy(j) 

&  /  (DYYS(J+1)  *  (DYYS(J)  +  DYYS(J+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CN)  -  CN) ) 

AS(I,J,K)  =  VISS1  +  (0.5  *  CS  *  DYY(J)  /  DYYS(J) 

&  +  0.125  *  (abs(cs)  +  cs)  *  dyy ( j )  *  dyy(j-l) 

&  /  (dyys(j-l)  *  dyys(j)) 

&  +  0.125  *  (abs(cs)  -  cs)  *  dyy(j-l)  *  dyy(}) 

&  /  (dyys(j)  *  (dyys(j)  +  dyys(j+l))) 

&  +  0.125  *  (abs(cn)  +  en)  *  dyy(j+l)  *  dyy(j) 

&  /  (DYYS(J)  *  (DYYS(J)  +  DYYS(J+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CS)  +  CS)) 

AF(I,J,K)  =  VISF1  +  (-  0.5  *  CF  *  DZZ(K)  /  DZZS(K+1) 

&  +  0.125  *  (abs(cf)  +  cf)  *  dzz(k)  *  dzz(k+l) 

&  /  (dzzs(k+l)  *  (dzzs(k)  +  dzzs(k+l))) 

&  +  0.125  *  (abs(cf)  -  cf)  *  dzz(k)  *  dzz(k+l) 

&  /  (dzzs(k+l)  *  dzzs(k+2)) 

&  +  0.125  *  (abs(cb)  -  cb)  *  dzz(k-l)  *  dzz(k) 

&  /  (DZZS(K+1)  *  (DZZS(K)  +  DZZS(K+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CF)  -  CF) ) 

AB(I,J,K)     =    VISB1    +     (0.5    *    CB    *    DZZ(K)     /    DZZS(K) 

&  +    0.125    *    (abs(cb)    +   cb)    *    dzz(k)    *    dzz(k-l) 

&  /    (dzzs(k-l)    *   dzzs(k) ) 

&  +    0.125    *    (abs(cb)    -   cb)    *    dzz(k-l)    *    dzz(k) 

&  /    (dzzs(k)    *    (dzzs(k)    +   dzzs(k+l))) 

&  +    0.125    *    (abs(cf)    +   cf)    *   dzz(k+l)    *    dzz(k) 

&  /     (DZZS(K)     *     (DZZS(K)     +    DZZS(K+1))))     *    QUICK 

&  +    UPWIND    *     (     0.5    *     (ABS(CB)     +    CB)) 

AEE(I,J,K)  =  -    0.0625    *    QUICK 

&  *  (ABS(CE)     -    CE)     *    DXX(I)**2 

&  /  (DXX(I+1)     *    DXXS(I+1) ) 

AWW(I,J,K)  =  -    0.0625    *    QUICK 

&  *  (ABS(CW)     +    CW)     *    DXX(I-1)**2 

&  /  (DXX(I-2)     *    DXXS(I-l) ) 

ANN(I,J,K)  =  -    0.125    *    QUICK 
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&  *  (ABS(CN)  -  CN)  *  DYY(J)  *  DYY(J+1) 

&  /  (DYYS(J+2)  *  (DYYS(J+1)  +  DYYS(J+2))) 

ASS  (I,  J.,  K)  =  -  0.125  *  QUICK 

&  *  (ABS(CS)  +  CS)  *  DYY(J)  *  DYY(J-l) 

&  /  (DYYS(J-l)  *  (DYYS(J-l)  +  DYYS(J))) 

AFF(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CF)  -  CF)  *  DZZ(K)  *  DZZ(K+1) 

&  /  (DZZS(K+2)  *  (DZZS(K+1)  +  DZZS(K+2))) 

ABB(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CB)  +  CB)  *  DZZ(K)  *  DZZ(K-l) 

&  /  (DZZS(K-l)  *  (DZZS(K-l)  +  DZZS(K))) 

AP(I,J,K)  =  AE(I,J,K)  +  AW(I,J,K)  +  AN(I,J,K)  +  AS(I,J,K) 

&  +  AF(I,J,K)  +  AB(I,J,K)  +  AWW(I,J,K) 

&  +  AEE(I,J,K)  +  ASS(I,J,K)  +  ANN (I, J, K) 

&  +  ABB (I, J, K)  +  AFF(I,J,K)  +  VOLDT 

SU(I,J,K)  =  VOLDT  *  UOD(I,J,K) 

&  +  DYZS  *  (P(I-1,J,K)  -  P(I,J,K)) 
100  continue 
C 
C      CALCULATION  OF  THE  SOURCE  TERM 

C 

DO  101  I=NUST+1,NI 
DO  101  J=2,NJ 

DO  101  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AWW(I,J,K)  *  UPD(I-2,J,K) 

101   CONTINUE 


C 


c 


c 


DO  102  I=NUST,NI-1 
DO  102  J=2,NJ 

DO  102  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AEE(I,J,K)  *  UPD(I+2,J,K) 

102  CONTINUE 

DO  103  I=NUST,NI 
DO  103  J=4,NJ 

DO  103  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ASS(I,J,K)  *  UPD(I,J-2,K) 

103  CONTINUE 

DO  104  I=NUST,NI 
DO  104  J=2,NJ-2 
DO  104  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ANN(I,J,K)  *  UPD(I,J+2,K) 

104  CONTINUE 

DO  105  I=NUST,NI 
DO  105  J=2,NJ 

DO  105  K=4,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ABB(I,J,K)  *  UPD(I,J,K-2) 

105  CONTINUE 


C 


C 


DO  106  I=NUST,NI 
DO  106  J=2,NJ 

DO  106  K=2,NK-2 
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SU(I,J,K)  -  SU(I,J,K) 
&  +  AFF(I,J,K)  *  UPD(I,J,K+2) 

106   CONTINUE 

c         BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 

c 

DO  121  K=2,NK 

DO  121  I=NUST,NI 

AP(I,2,K)  =  AP(I,2,K)  +  AS(I,2,K)  -  ASS(I,2,K) 
AS(I,2,K)  =  0. 

SU(I,3,K)  =  SU(I,3,K)  -  ASS(I,3,K)  *  UPD(I,2,K) 
SU(I,NJM1,K)  =  SU(I,NJM1,K) 
&  -  ANN(I,NJM1,K)  *  UPD(I,NJ,K) 

AP(I,NJ,K)  =  AP(I,NJ,K)  +  AN(I,NJ,K)  -  ANN(I,NJ,K) 
AN(I,NJ,K)  =  0. 

121  CONTINUE 
c 

c         LEFT  AND  RIGHT  WALL 

c 

DO  122  K=2,NK 

DO  122  J=2,NJ 

AW(NUST,J,K)  =  0. 
SU(NUST,J,K)  =  SU(NUST,J,K) 
&  -  AWW(NUST,J,K)  *  UPD(NUST,J,K) 

SU(NI,J,K)  =  SU(NI,J,K) 
&  -  AEE(NI,J,K)  *  UPD(NI,J,K) 

AE(NI,J,K)  =  0. 

122  CONTINUE 
c 

C         FRONT  AND  BACK  WALL 

c 

DO  12  3  J=2,NJ 

DO  123  I=NUST,NI 

AP(I,J,2)  =  AP(I,J,2)  +  AB(I,J,2)  -  ABB (I, J, 2) 
AB(I,J,2)  =  0. 

SU(I,J,3)  =  SU(I,J,3)  -  ABB (I, J, 3)  *  UPD(I,J,2) 
SU(I,J,NKM1)  =  SU(I,J,NKM1) 
&  -  AFF(I,J,NKM1)  *  UPD(I,J,NK) 

AP(I,J,NK)  =  AP(I,J,NK)  +  AF(I,J,NK)  -  AFF(I,J,NK) 
AF(I,J,NK)  =  0. 

123  CONTINUE 
c 

IF  (NCHP  .NE.  0)  THEN 
c 
c  U-VELOCITY  SET  TO  ZERO  IN  CHIP 

c 

DO  130  M=1,NJCHIP 

DO  130  N=1,NKCHIP 

DO  130  I=IBGN+1,IEND+1 

DO  130  J=JBGN(M) , JEND(M) 

DO  130  K=KBGN(N) ,KEND(N) 
AP(I,J,K)  =  small 
AW(I,J,K)  =  0. 
AE(I,J,K)  =  0. 
AS(I,J,K)  =  0. 
AN(I,J,K)  =  0. 
AB(I,J,K)  =  0. 
AF(I,J,K)  =  0. 
SU(I,J,K)  =  0. 
130      CONTINUE 
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C         CHIP  BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 

c 

DO  131  M=1,NJCHIP 

DO  131  N=1,NKCHIP 

DO  131  I=IBGN+1,IEND+1 

DO  131  K=KBGN(N) ,KEND(N) 
JJS  =  JEND(M) 
JJS1  =  JEND(M)  +  1 
JJS2  =  JEND(M)  +  2 
JJN  =  JBGN(M) 
JJN1  =  JBGN(M)  -  1 
JJN2  =  JBGN(M)  -  2 
AP(I,JJS1,K)  =  AP(I,JJS1,K) 
&  +  AS(I,JJS1,K)  -  ASS(I,JJS1,K) 

AS(I,JJS1,K)  =  0. 
SU(I,JJS2,K)  =  SU(I,JJS2,K) 
&  -  ASS(I,JJS2,K)  *  UPD(I,JJS,K) 

&  -  ASS(I,JJS2,K)  *  UPD(I,JJS1,K) 

SU(I,JJN2,K)  =  SU(I,JJN2,K) 
&  -  ANN(I, JJN2,K)  *  UPD(I,JJN,K) 

&  -  ANN(I, JJN2,K)  *  UPD(I,JJN1,K) 

AP(I,JJN1,K)  =  AP(I,JJN1,K) 
&  +  AN(I,JJN1,K)  -  ANN(I, JJN1,K) 

AN(I,JJN1,K)  =  0. 
131  CONTINUE 
c 

C         FRONT  AND  BACK  WALL 
c 

DO  132  M=1,NJCHIP 

DO  132  N=1,NKCHIP 

DO  132  I=IBGN+1,IEND+1 

DO  132  J=JBGN(M) , JEND(M) 
KKB  =  KEND(N) 
KKB1  =  KEND(N)  +  1 
KKB2  =  KEND(N)  +  2 
KKF  =  KBGN(N) 
KKF1  =  KBGN(N)  -  1 
KKF2  =  KBGN(N)  -  2 
AP(I,J,KKB1)  =  AP(I,J,KKB1) 
&  +  AB(I,J,KKB1)  -  ABB(I, J,KKB1) 

AB(I,J,KKB1)  =  0. 
SU(I,J,KKB2)  =  SU(I,J,KKB2) 
&  -  ABB(I,J,KKB2)  *  UPD(I,J,KKB) 

&  -  ABB(I,J,KKB2)  *  UPD ( I , J , KKB1 ) 

SU(I,J,KKF2)  =  SU(I,J,KKF2) 
&  -  AFF(I,J,KKF2)  *  UPD (I , J, KKF) 

&  -  AFF(I,J,KKF2)  *  UPD ( I , J , KKF1) 

AP(I,J,KKF1)  =  AP(I,J,KKF1) 
&  +  AF(I;J,KKF1)  -  AFF(I,J,KKF1) 

AF(I,J,KKF1)  =  0. 
132  CONTINUE 
C 

C  RIGHT  WALL 

c 

DO    133    M=1,NJCHIP 

DO    133    N=1,NKCHIP 

DO    133    J=JBGN(M) , JEND(M) 

DO    133    K=KBGN(N) ,KEND(N) 
112    =    IEND    +    2 
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AW(II2,J,K)  =  0. 

SU(II2,J,K)  =  SU(II2,J,K) 

&  -  AWW(II2,J,K)  *  UPD(IEND,J,K) 

&  -  AWW(II2,J,K)  *  UPD(II2,J,K) 
133     CONTINUE 
ENDIF 

SOLVE  THE  LINEARISED  U  MOMENTUM  EQUATIONS 

CALL  SIP  (NUST,2,2,NI,NJ;NK,U) 

CALCULATE  DU  NEEDED  FOR  PRESSURE  CORRECTION  AND  SOLVE 

DO  200  K=2,NK 

DO  200  J=2,NJ 

DO  200  I=NUST,NI 

SU(I,J,K)  =  DYY(J)  *  DZZ(K) 
200  CONTINUE 

c 

CALL  SIP  (NUST,2,2,NI,NJ,NK,DU) 

c 

END 
c* ********************************************************************* 

SUBROUTINE  CALV 

C 

C      THIS  SUBROUTINE  ASSEMBLES  THE  V-MOMENTUM  EQUATION  AND 

C      BOUNDARY  CONDITION  AND  SOLVES  FOR  THE  V-VELOCITY  (Y  DIRECTION) 

C 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ ra , pr , sorsum , dt ime , xper , rol 1 , dt_inv , xt  ime 

common/ limits/ni , nipl , niml ,  nj , n jpl , njml , nk, nkpl , nkml , 
&  nip2,njp2,nkp2, iter,nnmax 

common/fv_init/tod(12,16,16) , uod ( 12 , 16 , 16) ,vod(12, 16, 16) , 
&  wod(12, 16, 16) ,pod(12, 16, 16) 

common/fv_cur/t(12, 16,16),u(12,16,16),v(12,16,16), 
&  w(12,16,16) ,p(12,16,16) 

common/fv_int/tpd(12,16,16) , upd ( 12 , 16 , 16) , vpd ( 12 , 16 , 16) , 
&  wpd(12,16, 16) ,ppd(12, 16, 16) 

common/mscn/smp(12, 16, 16) ,resorm(93) , 
&  du(12, 16,  16)  ,dv(12, 16,16)  , 

&  dw(12, 16, 16) ,pp(12, 16, 16) 

common/coeff/ap(12, 16 , 16) , aw( 12 , 16, 16) , ae ( 12 , 16 , 16) , 
&  as ( 12, 16, 16 ), an (12, 16,16) ,ab( 12, 16,16) , 

&  af (12, 16, 16) ,su(12, 16, 16) 

COMMON/COEF2/AWW(12, 16, 16) ,AEE(12, 16, 16) , ASS (12, 16, 16) , 
&  ANN(12,16,16) f ABB (12, 16, 16) , AFF ( 12 , 16 , 16) 

common/condu/alpha(12,16,16) ,rho(12,16, 16) ,visco(12, 16, 16) 

common/chip/ ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) ,nchp, 
&  isub, ichip, jchip, kchip 

COMMON/SCHEME/QUICK, UPWIND 

common/dif f /ale , a Is , thot , tcool , tavg 

common/array/njehip, nkchip, ichoice 

common /tol/ small , eps, sormax 

common/ if ir st/ nust , nvst , nwst , npst 
c 

c         CALCULATE  COEFFICIENTS 
c 

do  100  k=2,nk 
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do  100  j=3,nj 

do  100  i=nvst,ni 

voldt  =  dxx(i)  *  dyys(j)  *  dzz(k)  /  dtime 

C        CALCULATE  INTERFACIAL  KINEMATIC  VISCOSITIES 
C  USING  THE  HARMONIC  MEAN  FORMULATION 

C  .  .  . 

visw  =  visco(i,j,k)  *  visco(i-l,],K) 

&  *  (dxx(i-l)  +  dxx(i))  /  (dxx(i-l) 

&  *  visco(i,j,k)  +  dxx(i)  *  visco ( l-l , j , k) ) 

vise  =  visco(i,j,k)  *  visco( i+l, j ,k) 
&  *  (dxx(i+l)  +  dxx(i))  /  (dxx(i+l) 

&  *  visco(i,j,k)  +  dxx(i)  *  visco( i+l, j , k) ) 

viss  =  visco(i, j-l,k) 

visn  =  visco(i,j,k) 

visb  =  visco(i,j,k)  *  visco( i , j ,k-l) 
&  *  (dzz(k-l)  +  dzz(k))  /  (dzz(k-l) 

&  *  visco(i,j ,k)  +  dzz(k)  *  visco( i , j , k-1) ) 

visf  =  visco(i,j ,k)  *  visco( i , j , k+i) 
&  *  (dzz(k+l)  +  dzz(k))  /  (dzz(k+l) 

&  *  visco(i,j,k)  +  dzz(k)  *  visco ( i , j , k+l) ) 

c 

dxys  =  dxx(i)  *  dyys(^) 
dyzs  =  dyys(j)  *  dzz(k) 
dzxs  =  dzz(k)  *  dxx(i) 
zxoyn  =  dzxs  /  dyy(j) 
zxoys  =  dzxs  /  dyy(j-l) 
xyozf  =  dxys  /  dzzs(k+l) 
xyozb  =  dxys  /  dzzs(k) 
yzoxe  =  dyzs  /  dxxs(i+l) 
yzoxw  =  dyzs  /  dxxs(i) 

c 

GN  =  V(I,J+1,K) 
GP  =  V(I,J,K) 
GS  =  V(I,J-1,K) 
GE  =  U(I+1,J,K) 
GSE  =  U(I+1,J-1,K) 
GW  =  U(I,J,K) 
GSW  =  U(I,J-1;K) 
GF  =  W(I,J,K+1) 
GSF  =  W(I,J-1,K+1) 
GB  =  W(I,J,K) 
GSB  =  W(I,J-1,K) 

c 

en  =  0.5  *  (gn  +  gp)  *  dzxs 

cs  =  0.5  *  (gp  +  gs)  *  dzxs 

ce  =  (ge  *  dyy(j-l)  +  gse  *  dyy(j)) 
&  /  (dyy(j-l)  +  dyy(j))  *  dyzs 

cw  =  (gw  *  dyy(j-l)  +  gsw  *  dyy(j)) 
&  /  (dyy(j-l)  +  c3yy(j))  *  dyzs 

cf  =  (gf  *  dyy(j-l)  +  gsf  *  dyy(j)) 
&  /  (dyy(j-i)  +  dyy(j))  *  dxys 

cb  =  (gb  *  dyy(j-l)  +  gsb  *  dyy(j)) 
&  /  (dyy(j-l)  +  dyy(j))  *  dxys 

c 

visel  =  yzoxe  *  vise 
viswl  =  yzoxw  *  visw 
visnl  =  zxoyn  *  visn 
vissl  =  zxoys  *  viss 
visfl  =  xyozf  *  visf 
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v  J.i3i_f  J.  —  ajt  ' 


AE(I  J  K)  =  VISE1  +  (-  0.5  *  CE  *  DXX(I)  /  DXXS(I+1) 

&  +  0.125  *  (abs(ce)  +  ce)  *  dxx(i)  *  dxx(i+l) 

&  /  (dxxs(i+l)  *  (dxxs(i)  +dxxs(i+l))) 

&  +  0.125  *  (abs(ce)  -  ce)  *  dxx(i)  *  dxx(i+l) 

&  /  (dxxs(i+l)  *  dxxs(i+2)) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (DXXS(I+1)  *  (DXXS(I)  +  DXXS(I+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CE)  -  CE) ) 

AW (I  J,K)  =  VISW1  +  (0.5  *  CW  *  DXX(I)  /  DXXS(I) 

&  +  0.125  *  (abs(cw)  +  cw)  *  dxx(i)  *  dxx(i-l) 

&  /  (dxxs(i-l)  *  dxxs(i)) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (dxxs(i)  *  (dxxs(i)  +dxxs(i+l))) 

&  +  0.125  *  (abs(ce)  +  ce)  *  dxx(i+l)  *  dxx(i) 

&  /  (DXXS(I)  *  (DXXS(I)  +  DXXS(I+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CW)  +  CW) ) 

AN(I  J,K)  =  VISN1  +  (-  0.5  *  CN  +  0.0625  *  (ABS(CN)  +  CN) 

&  *  dyy(j)  /  dyys(j)  +  0.125  *  (abs(cn)  -  en) 

&  *  dyy(j)  /  dyy(j+i)  +  0.0625  *  (abs(cs)  -  cs) 

&  *  DYY(J-1)**2  /  (DYY(J)  *  DYYS(J)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CN)  -  CN) ) 

AS(I  J  K)  =  VISS1  +  (0.5  *  CS  +  0.125  *  (ABS(CS)  +  CS) 

&  *  dyy(j-l)  /  dyy(j-2)  +  0.0625  *  (abs(cs)  -  cs) 

&  *  dyy(j-l)  /  dyys(j)  +  0.0625  *  (abs(cn)  +  en) 

&  *  DYY(J)**2  /  (DYY(J-l)  *  DYYS(J)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CS)  +  CS)) 

AF(I  J  K)  =  VISF1  +  (-  0.5  *  CF  *  DZZ(K)  /  DZZS(K+1) 

&  +  0.125  *  (abs(cf)  +  cf)  *  dzz(k)  *  dzz(k+l) 
&  /  (dzzs(k+l)  *  (dzzs(k)  +dzzs(k+l))) 

&  +  0.125  *  (abs(cf)  -  cf)  *  dzz(k)  *  dzz(k+l) 
&  /  (dzzs(k+i)  *  dzzs(k+2)) 

&  +  0.125  *  (abs(cb)  -  cb)  *  dzz(k-l)  *  dzz(k) 

&  /  (DZZS(K+1)  *  (DZZS(K)  +  DZZS(K+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CF)  -  CF) ) 

AB(I,J,K)  =  VISB1  +  (0.5  *  CB  *  DZZ(K)  /  DZZS(K) 
&  +  0.125  *  (abs(cb)  +  cb)  *  dzz(k)  *  dzz(k-l) 

&  /  (dzzs(k-l)  *  dzzs(k) ) 

&  +  0.125  *  (abs(cb)  -  cb)  *  dzz(k-l)  *  dzz(k) 

&  /  (dzzs(k)  *  (dzzs(k)  +  dzzs(k+i))) 

&  +  0.125  *  (abs(cf)  +  cf)  *  dzz(k+l)  *  dzz(k) 

&  /  (DZZS(K)  *  (DZZS(K)  +  DZZS(K+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CB)  +  CB) ) 

AEE(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CE)  -  CE)  *  DXX(I)  *  DXX(I+1) 

&  /  (DXXS(I+2)  *  (DXXS(I+1)  +  DXXS(I+2))) 

AWW(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CW)  +  CW)  *  DXX(I)  *  DXX(I-l) 

&  /  (DXXS(I)  *  (DXXS(I-l)  +  DXXS(I))) 

ANN(I,J,K)  =  -  0.0625  *  QUICK 

&  *  (ABS(CN)  -  CN)  *  DYY(J)**2 

&  /  (DYY(J+1)  *  DYYS(J+1)) 

ASS(I,J,K)  =  -  0.0625  *  QUICK 

&  *  (ABS(CS)  +  CS)  *  DYY(J-1)**2 

&  /  (DYY(J-2)  *  DYYS(J-l) ) 

AFF(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CF)  -  CF)  *  DZZ(K)  *  DZZ(K+1) 

&  /  (DZZS(K+2)  *  (DZZS(K+1)  +  DZZS(K+2))) 
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ABB(I,J,K)  =  -  0.125  *  QUICK 
&  *  (ABS(CB)  +  CB)  *  DZZ(K)  *  DZZ(K-l) 

&  /  (DZZS(K)  *  (DZZS(K-l)  +  DZZS(K))) 

AP(I,J,K)  -  AE(I,J,K)  +  AW(I,J,K)  +  AN(I,J,K)  +  AS(I,J,K) 

&  +  AF(I,J,K)  +  AB(I,J,K)  +  AWW(I,J,K) 

&  +  AEE(I,J,K)  +  ASS(I,J,K)  +  ANN(I,J,K) 

&  +  ABB (I, J, K)  +  AFF(I,J,K)  +  VOLDT 

SU(I,J,K)  =  VOLDT  *  VOD(I,J,K) 

&  +  DZXS  *  (P(I,J-1,K)  -  P(I,J,K)) 

&  +  RA  *  (T(I,J,K)  *  DYY(J-l)  +  T(I,J-1,K) 

&  *  DYY(J))  /  (DYY(J)  +  DYY(J-l)) 

&  *  DXX(I)  *  DYYS(J)  *  DZZ(K) 
100  CONTINUE 
C 
C      CALCULATION  OF  THE  SOURCE  TERM 

C 

DO  101  I=NVST+2,NI 
DO  101  J=3,NJ 

DO  101  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AWW(I,J,K)  *  VPD(I-2,J,K) 

101   CONTINUE 


C 


c 


DO  102  I=NVST,NI-2 
DO  102  J=3,NJ 

DO  102  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AEE(I,J,K)  *  VPD(I+2,J,K) 

102  CONTINUE 

DO  103  I=NVST,NI 
DO  103  J=4,NJ 

DO  103  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ASS(I,J,K)  *  VPD(I,J-2,K) 

103  CONTINUE 
C 

DO  104  I=NVST,NI 
DO  104  J=3,NJM1 
DO  104  K=2,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ANN(I,J,K)  *  VPD(I,J+2,K) 

104  CONTINUE 
C 

DO  105  I=NVST,NI 
DO  105  J=3,NJ 

DO  105  K=4,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ABB(I,J,K)  *  VPD(I,J,K-2) 

105  CONTINUE 
C 

DO  106  I=NVST,NI 
DO  106  J=3,NJ 

DO  106  K=2,NK-2 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AFF(I,J,K)  *  VPD(I,J,K+2) 

106  CONTINUE 
C 

C         BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM  WALL) 

C 
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DO    121    K=2,NK 

DO    121    I=NVST,NI 
AS(I,3,K)     =    0. 
SU(I,3,K)     =    SU(I,3,K) 
&  -    ASS(I,3,K)     *    VPD(I,3,K) 

SU(I,NJ,K)     =    SU(I,NJ,K) 
&  -    ANN(I,NJ,K)     *    VPD(I,NJ,K) 

AN(I,NJ,K)     =    0. 

121  CONTINUE 

C 

C  LEFT  AND  RIGHT  WALL 

C 

DO  122  K=2,NK 

DO  122  J=3,NJ 

AP(NVST,J,K)  =  AP(NVST,J,K)  +  AW(NVST,J,K) 
&  -  AWW(NVST, J,K) 

AW(NVST,J,K)  =  0. 

SU(NVST+1,J,K)  =  SU(NVST+1,J,K) 
&  -  AWW(NVST+1,J,K)  *  VPD(NVST,J,K) 

SU(NIM1,J,K)  =  SU(NIM1,J,K) 
&  -  AEE(NIM1,J,K)  *  VPD(NI,J,K) 

AP(NI,J,K)  =  AP(NI,J,K)  +  AE(NI,J,K) 

&  -  AEE(NI,J,K) 

AE(NI,J,K)  =  0. 

122  CONTINUE 
c 

c         FRONT  AND  BACK  WALL 

c 

DO  123  J=3,NJ 

DO  123  I=NVST,NI 

AP(I,J,2)  =  AP(I,J,2)  +  AB(I,J,2)  -  ABB (I, J, 2) 
SU(I,J,3)  =  SU(I,J,3) 
&  -  ABB (I, J, 3)  *  VPD(I,J,2) 

SU(I;J,NKM1)  =  SU(I,J,NKM1) 
&  -  AFF(I,J,NKM1)  *  VPD(I,J,NK) 

AB(I,J,2)  =  0. 

AP(IfJ,NK)  =  AP(I,J,NK)  +  AF(I,J,NK)  -  AFF(I,J,NK) 

AF(I,J,NK)  =  0. 

123  CONTINUE 
C 

IF  (NCHP  .NE.  0)  THEN 

c 

c  V-VELOCITY  SET  TO  ZERO  IN  CHIP 

c 

DO  130  M=1,NJCHIP 

DO  130  N=1,NKCHIP 

DO  130  I=IBGN,IEND 

DO  130  J=JBGN(M) ,JEND(M)+1 
DO  130  K=KBGN(N) ,KEND(N) 
AP(I,J,K)  =  small 
AW(I,J,K)  =  0.0 
AE(I,J,K)  =  0.0 
AS(I,J,K)  =  0.0 
AN(I,J,K)  =  0.0 
AB(I,J,K)  =  0.0 
AF(I,J;K)  =  0.0 
SU(I,J,K)  =  0.0 
130      CONTINUE 
c 
C         CHIP  BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 
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DO  131  M=1,NJCHIP 

DO  131  N=1,NKCHIP 

DO  131  I=IBGN,IEND 

DO  131  K=KBGN(N) ,KEND(N) 
JJS  =  JEND(M) 
JJS2  =  JEND(M) 
JJN1  =  JBGN(M) 
JJP1  =  JBGN(M) 
AS(I,JJS2,K)  = 
SU(I,JJS2,K)  = 


& 
& 

& 
& 

131  CONTINUE 


SU(I,JJN1,K)  = 


AN(I,JJN1,K)  = 


+  2 

-  1 

+  1 

0. 

SU(I, JJS2,K) 

ASS(I, JJS2,K) 

ASS(I, JJS2,K) 

SU(I, JJN1,K) 

ANN(I,JJN1,K) 

ANN(I, JJN1,K) 

0. 


*  VPD(I,JJS,K) 

*  VPD(I, JJS2,K) 

*  VPD(I, JJP1,K) 

*  VPD(I,JJN1,K) 


c 
C 

c 


FRONT  AND  BACK  WALL 

DO  132  M=1,NJCHIP 

DO  132  N=1,NKCHIP 

DO  132  I=IBGN,IEND 

DO  132  J=JBGN(M) , JEND(M)+1 
KKB  =  KEND(N) 


KKB1  =  KEND(N) 

+  1 

KKB2  =  KEND(N) 

+  2 

KKF  =  KBGN(N) 

KKF1  =  KBGN(N) 

-  1 

KKF2  =  KBGN(N) 

-  2 

AP(I,J,KKB1)  = 

AP(I, J,KKB1) 

& 

+ 

AB(I, J,KKB1) 

AB(I,J,KKB1)  = 

0. 

SU(I,J,KKB2)  = 

SU(I, J,KKB2) 

& 

- 

ABB(I, J,KKB2) 

& 

- 

ABB(I, J,KKB2) 

SU(I,J,KKF2)  = 

SU(I, J,KKF2) 

& 

- 

AFF(I, J,KKF2) 

& 

- 

AFF(I, J,KKF2) 

AP(I,J,KKF1)  = 

AP(I, J,KKF1) 

& 

+ 

AF(I, J,KKF1) 

AF(I,J,KKF1)  = 

0. 

132  CONTINUE 

c 
C 

RIGHT  WALL 

c 

-  ABB(I, J,KKB1) 


*  VPD(I,J,KKB) 

*  VPD(I, J,KKB1) 

*  VPD(I,J,KKF) 

*  VPD(I, J,KKF1) 

-  AFF(I,J,KKF1) 


DO    133    M=1,NJCHIP 

DO    133    N=1,NKCHIP 

DO    133    J=JBGN(M) , JEND(M) 

DO    133    K=KBGN(N) ,KEND(N) 

II  =    IEND 

III  =  IEND  +  1 
112  =  IEND  +  2 
AP(II1,J,K)     =    AP(II1,J,K) 

+ 
SU(II2,J,K)     = 


AW(II1, J,K) 


AW(II1, J,K) 
SU(II2, J,K) 
AWW(II2, J,K) 
AWW(II2, J,K) 
0. 


-    AWW(II1,J,K) 


VPD(II, J,K) 
VPD(II1, J,K) 
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200 
c 

c 

c**** 

c 

C 

c 
c 


c 
c 

c 


CONTINUE 
ENDIF 

SOLVE  THE  LINEARISED  V  MOMENTUM  EQUATION 
CALL  SIP  (NVST,3,2,NI,NJ,NK,V) 

SET  UP  SU  FOR  CALCULATING  DV  AND  SOLVE 

DO  200  K=2,NK 

DO  200  J=3,NJ 

DO  200  I=NVST,NI 

SU(I,J,K)  =  DZZ(K)  *  DXX(I) 
CONTINUE 

CALL  SIP  (NVST,3,2,NI,NJ,NK,DV) 

END 
******************************************************************* 

SUBROUTINE  CALW 

THIS  SUBROUTINE  ASSEMBLES  THE  W-MOMENTUM  EQUATION  AND 
BOUNDARY  CONDITION  AND  SOLVES  FOR  THE  W-VELOCITY  (Z  DIRECTION) 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ra , pr , sorsum , dtime , xper , roll , dt_inv , xt line 

common/ limits/ni , nipl , niml , nj , njpl , n jml , nk, nkpl , nkml , 
&  nip2,njp2,nkp2, iter,nnmax 

common/fv_init/tod(12,16,16) , uod ( 12 , 16, 16) , vod ( 12 , 16 , 16) , 
&  wod(12,16,16) ,pod(12,16,16) 

common/fv_cur/t(12,16,16) ,u(12,16,16),v(12,16,16), 
&  w(12,16,16) ,p(12, 16,16) 

common/ fv_int/tpd (12,16,16), upd (12,16 
&  wpd(12,16,16) ,ppd(12,16 

common/mscn/smp(12, 16, 16) ,resorm(93) , 
&  du(12,16,16) ,dv(12,16,16) , 

&  dw(12, 16,16) ,pp(12,16,16) 

common/coeff/ap(12,16,16) , aw( 12 , 16, 16) , ae ( 12 , 16 , 16) , 
&  as (12, 16, 16) , an (12, 16,16) ,ab( 12, 16, 16) , 

&  af (12, 16, 16) ,su(12, 16, 16) 

COMMON/COEF2/AWW(12, 16, 16) ,AEE(12, 16, 16) , ASS (12, 16, 16) , 
&  ANN(12,16,16) ,ABB(12,16,16) , AFF ( 12 , 16 , 16) 

common/condu/alpha (12,16,16), rho (12,16,16) , visco (12,16,16) 

common/chip/ ibgn,iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) ,nchp, 
&  isub, ichip, jchip,kchip 

COMMON/SCHEME/QUICK, UPWIND 

common/diff /ale, als, thot , tcool, tavg 

common/array/njchip,nkchip, ichoice 

common /tol/ small, eps, sormax 

common/ if irst/nust , nvst , nwst , npst 


16) , vpd(12, 16 
16) 


16) 


CALCULATE  COEFFICIENTS 


do 


100  k=3,nk 
do  100  j=2,nj 

do  100  i=nwst,ni 
voldt  =  dxx(i) 


*  dyy(j)  *  dzzs(k)  /  dtime 
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C        CALCULATE  INTERFACIAL  KINEMATIC  VISCOSlTlfcb 
C  USING  THE  HARMONIC  MEAN  FORMULATION 


C 


visw  =  visco(i, j ,k)  *  visco( i-1 , j , k) 
&  *  (dxx(i-l)  +  dxx(i))  /  (dxx(i-l) 

&  *  visco(i,j,k)  +  dxx(i)  *  visco( i-1, j , k) ) 

vise  =  visco(i, j ,k)  *  visco( i+1 , j , k) 
&  *  (dxx(i+l)  +  dxx(i))  /  (dxx(i+l) 

&  *  visco(i,j,k)  +  dxx(i)  *  visco( i+1 , j , k) ) 

viss  =  visco(i,j,k)  *  visco(i, j-1 , k) 
&  *  (dyy(j-l)  +dyy(j))  /  (dyy(j-l) 

&  *  visco(i,j,k)  +  dyy(j)  *  visco( i , j-1 , k) ) 

visn  =  visco(i,j,k)  *  visco( i , j+1 , k) 
&  *  (dyy(j+l)  +dyy(j))  /  (dyy(j+l) 

&  *  visco(i,j,k)  +  dyy(j)  *  visco( i, j+1 , k) ) 

visb  =  visco(i, j ,k-l) 

visf  =  visco(i,j,k) 

dzxs  =  dzzs(k)  *  dxx(i) 
dxys  =  dxx(i)  *  dyy ( j ) 
dyzs  =  dyy(j)  *  dzzs(k) 
zxoyn  =  dzxs  /  dyys(j+l) 
zxoys  =  dzxs  /  dyys(j) 
xyozf  =  dxys  /  dzz(k) 
xyozb  =  dxys  /  dzz(k-l) 
yzoxe  =  dyzs  /  dxxs(i+l) 
yzoxw  =  dyzs  /  dxxs(i) 

GN  =  V(I,J+1,K) 
GNB  =  V(I,J+1,K-1) 
GS  =  V(I,J,K) 
GSB  =  V(I,J,K-1) 
GE  =  U(I+1,J,K) 
GEB  =  U(I+1,J,K-1) 
GW  =  U(I,J,K) 
GWB  =  U(I,J,K-1) 
GF  =  W(I,J,K+1) 
GP  =  W(I,J,K) 
GB  =  W(I,J,K-1) 

en  =  (gn  *  dzz(k-l)  +  gnb  *  dzz(k)) 

&  /  (dzz(k-l)  +  dzz(k))  *  dzxs 

cs  =  (gs  *  dzz(k-l)  +  gsb  *  dzz(k)) 

&  /  (dzz(k-l)  +  dzz(k))  *  dzxs 

ce  =  (ge  *  dzz(k-l)  +  geb  *  dzz(k)) 

&  /  (dzz(k-l)  +  dzz(k))  *  dyzs 

cw  =  (gw  *  dzz(k-l)  +  gwb  *  dzz(k)) 

&  /  (dzz(k-l)  +  dzz(k))  *  dyzs 

cf  =  0.5  *  (gf  +  gp)  *  dxys 

cb  =  0.5  *  (gb  +  gp)  *  dxys 

visel  =  yzoxe  *  vise 

viswl  =  yzoxw  *  visw 

visnl  =  zxoyn  *  visn 

vissl  =  zxoys  *  viss 

visfl  =  xyozf  *  visf 

visbl  =  xyozb  *  visb 

AE(I,J,K)     -    VISE1    +     (-    0.5    *    CE    *    DXX(I)     /    DXXS(I+1) 
&  +    0.125    *    (abs(ce)    +   ce)    *    dxx(i)    *   dxx(i+l) 
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&  /  (dxxs(i+l)  *  (dxxs(i)  +  dxxsii+ijjj 

&  +  0.125  *  (abs(ce)  -  ce)  *  dxx(i)  *  dxx(i+l) 

&  /  (dxxs(i+l)  *  dxxs(i+2)) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (DXXS(I+1)  *  (DXXS(I)  +  DXXS(I+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CE)  -  CE) ) 

AW(I,J,K)  =  VISW1  +  (0.5  *  CW  *  DXX(I)  /  DXXS(I) 

&  +  0.125  *  (abs(cw)  +  cw)  *  dxx(i)  *  dxx(i-l) 

&  /  (dxxs(i-l)  *  dxxs(i)) 

&  +  0.125  *  (abs(cw)  -  cw)  *  dxx(i-l)  *  dxx(i) 

&  /  (dxxs(i)  *  (dxxs(i)  +  dxxs(i+l))) 

&  +  0.125  *  (abs(ce)  +  ce)  *  dxx(i+l)  *  dxx(i) 

&  /  (DXXS(I)  *  (DXXS(I)  +  DXXS(I+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CW)  +  CW) ) 

AN(I,J,K)  =  VISN1  +  (-  0.5  *  CN  *  DYY(J)  /  DYYS(J+1) 

&  +  0.125  *  (abs(cn)  +  en)  *  dyy(j)  *  dyy(j+l) 

&  /  (dyys(j+l)  *  (dyys(j)  +  dyys(j+l))) 

&  +  0.125  *  (abs(cn)  -  en)  *  dyy ( j )  *  dyy(j+l) 

&  /  (dyys(j+l)  *  dyys(j+2)) 

&  +  0.125  *  (abs(cs)  -  cs)  *  dyy(j-l)  *  dyy(}) 

&  /  (DYYS(J+1)  *  (DYYS(J)  +  DYYS(J+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CN)  -  CN)) 

AS(I,J,K)  =  VISS1  +  (0.5  *  CS  *  DYY(J)  /  DYYS(J) 

&  +  0.125  *  (abs(cs)  +  cs)  *  dyy  ( j )  *  dyy(j-l) 

&  /  (dyys(j-l)  *  dyys(j) ) 

&  +  0.125  *  (abs(cs)  -  cs)  *  dyy(j-l)  *  dyy (j ) 

&  /  (dyys(j)  *  (dyys(j)  +  dyys(j+l))) 

&  +  0.125  *  (abs(cn)  +  en)  *  dyy(j+l)  *  dyy ( j ) 

&  /  (DYYS(J)  *  (DYYS(J)  +  DYYS(J+1))))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CS)  +  CS)) 

AF(I,J,K)  =  VISF1  +  (-  0.5  *  CF  +  0.0625  *  (ABS(CF)  +  CF) 

&  *  dzz(k)  /  dzzs(k)  +  0.125  *  (abs(cf)  -  cf) 

&  *  dzz(k)  /  dzz(k+l)  +  0.0625  *  (abs(cb)  -  cb) 

&  *  DZZ(K-1)**2  /  (DZZ(K)  *  DZZS(K)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CF)  -  CF) ) 

AB(I,J,K)  =  VISB1  +  (0.5  *  CB  +  0.125  *  (ABS(CB)  +  CB) 

&  *  dzz(k-l)  /  dzz(k-2)  +  0.0625  *  (abs(cb)  -  cb) 

&  *  dzz(k-l)  /  dzzs(k)  +  0.0625  *  (abs(cf)  +  cf) 

&  *  DZZ(K)**2  /  (DZZ(K-l)  *  DZZS(K)))  *  QUICK 

&  +  UPWIND  *  (  0.5  *  (ABS(CB)  +  CB) ) 

AEE(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CE)  -  CE)  *  DXX(I)  *  DXX(I+1) 

&  /  (DXXS(I+2)  *  (DXXS(I+1)  +  DXXS(I+2))) 

AWW(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CW)  +  CW)  *  DXX(I)  *  DXX(I-l) 

&  /  (DXXS(I-l)  *  (DXXS(I-l)  +  DXXS(I))) 

ANN(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CN)  -  CN)  *  DYY(J)  *  DYY(J+1) 

&  /  (DYYS(J+2)  *  (DYYS(J+1)  +  DYYS(J+2))) 

ASS(I,J,K)  =  -  0.125  *  QUICK 

&  *  (ABS(CS)  +  CS)  *  DYY(J)  *  DYY(J-l) 

&  /  (DYYS(J-l)  *  (DYYS(J-l)  +DYYS(J))) 

AFF(I,J,K)  =  -  0.0625  *  QUICK 

&  *  (ABS(CF)  -  CF)  *  DZZ(K)**2 

&  /  (DZZ(K+1)  *  DZZS(K+1)) 

ABB(I,J,K)  =  -  0.0625  *  QUICK 

&  *  (ABS(CB)  +  CB)  *  DZZ(K-1)**2 

&  /  (DZZ(K-2)  *  DZZS(K-l)) 

AP(I,J,K)  =  AE(I,J,K)  +  AW(I,J,K)  +  AN(I,J,K)  +  AS(I,J,K) 
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&  +  AF(I,J,K)  +  AB(I,J,K)  +  AWW(I,J,K) 

&  +  AEE(I,J,K)  +  ASS(I,J,K)  +  ANN(I,J,K) 

&  +  ABB (I, J, K)  +  AFF(I,J,K)  +  VOLDT 

SU(I,J,K)  =  VOLDT  *  WOD(I,J,K) 
&  +  DXYS  *  (P(I,J,K-1)  -  P(I,J,K)) 

100  continue 
C 

C      CALCULATION  OF  THE  SOURCE  TERM 
C 

DO  101  I=NWST+2,NI 
DO  101  J=2,NJ 

DO  101  K=3,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AWW(I,J,K)  *  WPD(I-2,J,K) 

101  CONTINUE 
C 

DO  102  I=NWST,NI-2 
DO  102  J=2,NJ 

DO  102  K=3,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AEE(I,J,K)  *  WPD(I+2,J,K) 

102  CONTINUE 
C 

DO  103  I=NWST,NI 
DO  103  J=4,NJ 

DO  103  K=3,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ASS(I,J,K)  *  WPD(I,J-2,K) 

103  CONTINUE 
C 

DO  104  I=NWST,NI 
DO  104  J=2,NJ-2 
DO  104  K=3,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ANN(I,J,K)  *  WPD(I,J+2,K) 

104  CONTINUE 
C 

DO  105  I=NWST,NI 
DO  105  J=2,NJ 

DO  105  K=4,NK 

SU(I,J,K)  =  SU(I,J,K) 
&  +  ABB(I,J,K)  *  WPD(I,J,K-2) 

105  CONTINUE 
C 

DO  106  I=NWST,NI 
DO  106  J=2,NJ 

DO  106  K=3,NKM1 

SU(I,J,K)  =  SU(I,J,K) 
&  +  AFF(I,J,K)  *  WPD(I,J,K+2) 

106  CONTINUE 
c 

C  BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 

C 

DO  121  K=3,NK 

DO  121  I=NWST,NI 

AP(I,2,K)  =  AP(I,2,K)  +  AS(I,2,K)  -  ASS(I,2,K) 
SU(I,3,K)  =  SU(I,3,K)  -  ASS(I,3,K)  *  WPD(I,2,K) 
SU(I,NJM1,K)  =  SU(I,NJM1,K) 
&  -  ANN(I,NJM1,K)  *  WPD(I,NJ,K) 

AS(I,2,K)  =  0. 
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AP(I,NJ,K)     =    AP(I,NJ,K)     +    AN(I,NJ,K)     -    ANN(I,NJ,K) 
AN(I,NJ,K)     =    0. 

121  CONTINUE 
c 

c  LEFT  AND  RIGHT  WALL 

c 

DO  122  K=3,NK 

DO  122  J=2,NJ 

AP(NWST,J,K)  =  AP(NWST,J,K)  +  AW(NWST,J,K) 
&  -  AWW(NWST, J,K) 

AW(NWST,J,K)  =  0. 

SU(NWST+1,J,K)  =  SU(NWST+1, J,K) 
&  -  AWW(NWST+1, J,K)  *  WPD(NWST, J,K) 

SU(NIM1,J,K)  =  SU(NIM1,J,K) 
&  -  AEE(NIM1,J,K)  *  WPD(NI,J,K) 

AP(NI,J,K)  =  AP(NI,J,K)  +  AE(NI,J,K)  -  AEE(NI,J,K) 
AE(NI,J,K)  =  0. 

122  CONTINUE 
c 

c  FRONT  AND  BACK  WALL 

c 

DO  12  3  J=2,NJ 

DO  123  I=NWST,NI 
AB(I, J, 3)  =  0. 

SU(I,J,3)  =  SU(I,J,3)  -  ABB(I,J,3)  *  WPD(I,J,3) 
SU(I,J,NK)  =  SU(I,J,NK) 
&  -  AFF(I,J,NK)  *  WPD(I,J,NK) 

AF(I,J,NK)  =  0. 
123   CONTINUE 
C 

if  (nchp  .ne.  0)  then 
c 

c  W-VELOCITY  SET  TO  ZERO  IN  CHIP 

c 

DO  130  M=1,NJCHIP 

DO  130  N=1,NKCHIP 

DO  130  I=IBGN,IEND 

DO  130  J=JBGN(M) , JEND(M) 

DO  130  K=KBGN(N) ,KEHD(N)+1 
AP(I,J,K)  =  small 
AW(I,J,K)  =  0.0 
AE(I,J,K)  =  0.0 
AS(I,J,K)  =  0.0 
AN(I,J,K)  =  0.0 
AB(I,J,K)  =  0.0 
AF(I,J,K)  =  0.0 
SU(I,J;K)  =  0.0 
130      CONTINUE 
c 

C         CHIP  BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 
c 

DO  131  M=1,NJCHIP 

DO  131  N=1,NKCHIP 

DO  131  I=IBGN,IEND 

DO  131  K=KBGN(N) ,KEND(N)+1 
JJS  =  JEND(M) 
JJS1  =  JEND(M)  +  1 
JJS2  =  JEND(M)  +  2 
JJN  =  JBGN(M) 
JJN1  =  JBGN(M)  -  1 
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JJN2    =    JBGN(M)     -    2 

AP(I,JJS1,K)     =    AP(I,JJS1,K) 
&  +    AS(I,JJS1,K)     -    ASS(I, JJS1,K) 

AS(I,JJS1,K)     =    0. 

SU(I,JJS2,K)     =    SU(I,JJS2,K) 
&  -    ASS(I,JJS2,K)     *    WPD(I,JJS,K) 

&  -    ASS(I,JJS2,K)     *    WPD(I,JJS1,K) 

SU(I,JJN2,K)     =    SU(I,JJN2,K) 
&  -    ANN (I, JJN2,K)     *    WPD(I,JJN,K) 

&  -    ANN(I,JJN2,K)     *    WPD(I,JJN1,K) 

AP(I,JJN1,K)     =    AP(I,JJN1,K) 
&  +    AN(I,JJN1,K)     -    ANN(I, JJN1,K) 

AN(I,JJN1,K)     =    0. 

131  CONTINUE 
C 

C         FRONT  AND  BACK  WALL 
c 

DO  132  M=1,NJCHIP 

DO  132  N=1,NKCHIP 

DO  132  I=IBGN,IEND 

DO  132  J=JBGN(M) , JEND(M) 
KKB  =  KEND(N) 
KKB2  =  KEND(N)  +  2 
KKP1  =  KBGN(N)  +  1 
KKF1  =  KBGN(N)  -  1 
AB(I,J,KKB2)  =  0. 
SU(I,J,KKB2)  =  SU(I,J,KKB2) 
&  -  ABB(I,J,KKB2)  *  WPD(I,J,KKB) 

&  -  ABB(I,J,KKB2)  *  WPD ( I , J , KKB2 ) 

SU(I,J,KKF1)  =  SU(I,J,KKF1) 
&  -  AFF(I,J,KKF1)  *  WPD(I,J,KKP1) 

&  -  AFF(I,J,KKF1)  *  WPD(I,J,KKF1) 

AF(I,J,KKF1)  =  0. 

132  CONTINUE 
c 

C  RIGHT  WALL 

c 

DO  133  M=1,NJCHIP 

DO  133  N=1,NKCHIP 

DO  133  J=JBGN(M) , JEND(M) 

DO  133  K=KBGN(N) ,KEND(N)+1 

II  =  IEND 

III  =  IEND  +  1 
112  =  IEND  +  2 
AP(II1,J,K)  =  AP(II1,J,K) 

&  +  AW(II1,J,K)  -  AWW(II1,J,K) 

SU(II2,J,K)  =  SU(II2,J,K) 
&  -  AWW(II2,J,K)  *  WPDCII^^) 

&  -  AWW(II2,J,K)  *  WPD(II1,J,K) 

AW(II1,J,K)  =  0. 

133  CONTINUE 
ENDIF 

c 

c      SOLVE  THE  LINEARISED  W  MOMENTUM  EQUATIONS 

c 

CALL  SIP  (NWST,2,3,NI,NJ,NK,W) 
C 

c      SET  UP  SU  FOR  CALCULATING  DW  AND  SOLVE 
c 

DO  200  K=3,NK 
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DO  200  J=2,NJ 

DO  200  I=NWST,NI 

su(i,j,k)  =  dxx(i)  *  dyy ( j ) 
200  CONTINUE 

c 

CALL  SIP  (NWST,2,3,NI,NJ,NK,DW) 

c 

END 
c*  ********************************************************************** 

SUBROUTINE  CALP 

C 

C      THIS  SUBROUTINE  ASSEMBLES  THE  PRESSURE  CORRECTION  EQUATION. 

C      THE  PRESSURE  CORRECTION  IS  SET  TO  ZERO  AT  ALL  THE  BOUNDARIES 

C 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ra , pr , sorsum, dtime , xper , roll , dt_inv, xtime 

common / to 1/ small , eps, sormax 

cornmon/dims/h,wth,bth,hchip,wchip,bchip,bsub 

common/ limits/ni , nipl , niml , nj , njpl , njml, nk, nkpl , nkml , 
&  nip2 ,njp2 , nkp2 , iter , nnmax 

common/ fv_init/tod( 12, 16, 16) ,uod(12, 16,16) , vod(12, 16, 16) , 
&  wod(12, 16, 16) ,pod(12, 16, 16) 

common/fv_cur/t(12,16,16) ,u(12,16,16),v(12,16,16), 
&  w(12, 16, 16) ,p(12, 16, 16) 

common/ fv_int/tpd( 12, 16, 16) ,upd(12, 16, 16) ,vpd(12, 16, 16) , 
&  wpd(12, 16, 16) ,ppd(12, 16, 16) 

common/mscn/smp(12, 16, 16) ,resorm(93) , 
&  du(12, 16,16) ,dv(12, 16,16) , 

&  dw(12, 16,16) ,pp(12,16, 16) 

common/coeff/ap(12,16, 16), aw (12, 16, 16), ae (12, 16, 16), 
&  as(12,16,16) , an (12, 16,16 ) ,ab( 12, 16,16 ) , 

&  af (12, 16,  16)  ,su(12,16,16) 

COMMON/COEF2/AWW(12, 16, 16) ,AEE(12, 16, 16) , ASS (12, 16, 16) , 
&  ANN(12,16,16) , ABB (12 ,16, 16) , AFF ( 12 , 16 , 16 ) 

common/mean/t_mean(12, 16, 16) ,u_mean(12, 16, 16) , 
&  v_mean(12, 16, 16) ,w_mean(12, 16, 16) , 

&  p_mean(12 , 16, 16) 

COMMON/SCHEME/QUICK, UPWIND 

common/array/njchip, nkchip, ichoice 

common/ if irst/nust , nvst , nwst, npst 


c 

c 

SET  UP  COEFFICIENTS 

c 

do  100  k=2,nk 

do  100  j=2,nj 

do  100  i=npst,ni 

dxys  =  dxx(i)  *  dyy(j) 

dyzs  =  dyy(j)  *  dzz(k) 

dzxs  =  dzz(k)  *  dxx(i) 

an(i,j,k)  =  dzxs  *  dv(i,j+l,k) 

as(i,j,k)  =  dzxs  *  dv(i,j,k) 

ae(i,j,k)  =  dyzs  *  du(i+l,j,k) 

aw(i,j,k)  =  dyzs  *  du(i,j,k) 

af(i,j,k)  =  dxys  *  dw(i,j,k+l) 

ab(i,j,k)  =  dxys  *  dw(i,j,k) 

en  =  v(i,j+l,k)  *  dzxs 

cs  =  v(i,j,k)  *  dzxs 

ce  =  u(i+l, j ,k)  *  dyzs 
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cw  =  u(i, j,k)  *  dyzs 

cf  =  w(i,j,k+l)  *  dxys 

cb  =  w(i,j,k)  *  dxys 

smp(i,j,k)  =  -  ce  +  cw  -  en  +  cs  -  cf  +  cb 

su(i, j,k)  =  smp(i, j ,k) 

100  continue 
c 

c  BOUNDARY  CONDITIONS  (TOP  AND  BOTTOM) 

c 

do  101  k=2,nk 

do  101  i=npst,ni 
as(i,2,k)  =  0. 
an(i,nj ,k)  =  0. 

101  continue 
c 

c  LEFT  AND  RIGHT  WALL 

c 

do  102  k=2,nk 
do  102  j=2,nj 

aw(npst,j,k)  =  0. 
ae (ni, j  ,k)  =0. 

102  continue 
c 

C  FRONT  AND  BACK  WALL 

c 

do  103  i=npst,ni 
do  103  j=2,nj 

ab(i,j,2)  =  0. 
af(i,j,nk)  =  0. 
103   continue 
c 

c        CALCULATE  AP 
c 

do  200  j=2,nj 

do  200  i=npst,ni 
do  200  k=2,nk 

ap(i,j,k)  =  an(i,j,k)  +  as(i,j,k)  +  ae(if j,k) 
&  +  aw(i,j,k)  +  af(i,j,k)  +  ab(i,j,k) 

200  continue 
c 

c        SOLVE  THE  PRESSURE  CORRECTION  EQUATION 

c 

call  sip  (npst,2,2,ni,nj ,nk,pp) 
c 

c        U  VELOCITY  CORRECTION 
c 

do  201  i=nust,ni 
do  201  j=2,nj 

do  201  k=2,nk 
U(I,J,K)  =  U(I,J,K)  +  DU(I,J,K)  *  (PP(I-1,J,K)  -  PP(I,J,K)) 

201  continue 
c 

c        V  VELOCITY  CORRECTION 
c 

do  202  j=3,nj 

do  202  k=2,nk 

do  202  i=nvst,ni 
V(I,J,K)  -  V(I,J,K)  +  DV(I,J,K)  *  (PP(I,J-1,K)  -  PP(I,J,K)) 

202  continue 
c 
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c       W  VELOCITY  CORRECTION 

C 

do  203  k=3,nk 

do  203  i=nwst,ni 
do  203  j=2,nj 
W(I,J,K)  =  W(I,J,K)  +  DW(I,J,K)  *  (PP(I,J,K-1)  -  PP(I,J,K)) 

203   continue 
c 
c       PRESSURE  CORRECTION 

c 

do  204  j=2,nj 

do  204  i=npst,ni 
do  204  k=2,nk 

p(i,j,k)  -  p(i,j,M  +  pp(i,j,k) 
PP(I,J,K)  =  0. 
204  continue 

RECALCULATE  MASS  FLUX  ERRORS  AFTER  U , V, W , P , CORRECTIONS 


c 
c 


sorsum  =  0 . 
resorm( iter)  =  0. 
do  205  j=2,nj 

do  205  i=npst,ni 
do  205  k=2,nk 

dxys  =  dxx(i)  *  dyy(j) 

dyzs  =  dyy(j)  *  dzz(k) 

dzxs  =  dzz(k)  *  dxx(i) 

en  =  v(i, j+l,k)  *  dzxs 

cs  =  v(i, j,k)  *  dzxs 

ce  =  u(i+l, j,k)  *  dyzs 

cw  =  u(i, j,k)  *  dyzs 

cf  =  w(i, j,k+l)  *  dxys 

cb  =  w(i, j ,k)  *  dxys 

smp(i,j,k)  =  -  ce  +  cw  -  en  +  cs  -  cf  +  cb 

sorsum  =  sorsum  +  smp(i,j,k) 

resorm(iter)  =  resorm(iter)  +  abs (smp ( i , j , k) ) 
205  continue 
end 
c  ********************************************************************* 

subroutine  nu 
C 
C     THIS  SUBROUTINES  DETERMINES  THE  OVERALL  HEAT  BALANCE  FOR  THE 

C     ENTIRE  ENCLOSURE 


C 


implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz (40) ,x(40) ,y (40) , z (40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common /parm/ra , pr , sor sum, dtime, xper , roll , dt_inv, xtime 

common / to 1/ small , eps, sormax 

common /dims/h, wth, bth, hchip, wchip, bchip, bsub 

common/ count /nt,nmax, imax, itmax, krun, nprint 

common/mscn/smp(12/ 16, 16) ,resorm(93) , 
&  du(12, 16,16) ,dv(12,16, 16) , 

&  dw(12,16,16) ,pp(12, 16,16) 

common/ limits/ni , nipl , niml , nj , njpl, njml , nk, nkpl , nkml , 
&  nip2 ,njp2 ,nkp2 , iter,nnmax 

common/fv_init/tod(12,16,16) ,uod ( 12 , 16 , 16) , vod ( 12 , 16 , 16) , 
&  wod(12, 16, 16) ,pod(12, 16, 16) 

common/ fv_cur/t( 12, 16,16),u(12,16,16),v(12,16,16), 
&  w(12, 16,16) ,p(12,16, 16) 
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common/ conau/ a ipna  (i^(ib,ioj  ,  rno^  i^id,idj  ,  v  x^^^\  ±*. ,  x<-> ,  j.^  , 
common/chip/ ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) , nchp, 
&  isub, ichip, jchip,kchip 

common/diff/alc,als,thot,tcool,tavg 
common/array/njchip, nkchip, ichoice 

c 

RNUC  =  0. 
RNUH  =  0. 
QCHIP  =  0. 
guns  =  0 . 
QCRT  =  0. 
do  10  i=l,ni 
do  10  j=l;nj 

do  10  k=l,nk 

quns  =  quns  +  (t(i,j,k)  -  tod(i, j ,k) )  *  dxx(i) 
&  *  dyy(j)  *  dzz(k)  /  dtime 

QCRT  =  QCRT  +  T(I,J,K)  *  SMP(I,J,K) 
10   continue 

al  = ' 1.  /  dyys(2) 
a2  =  1.  /  dyys(njpl) 
DO  20  1=2, NI 
DO  2  0  K=2,NK 

DTDYH  =  al  *  (T(I,2,K)  -  T(I,1,K))  *  ALPHA(I,2,K) 
DTDYC  =  a2  *  (T(I;NJP1,K)  -  T(I,NJ,K))  *  ALPHA ( I , NJ , K) 
RNUC  =  RNUC  -  DTDYC  *  dzz(k)  *  dxx(i) 
RNUH  =  RNUH  -  DTDYH  *  dzz(k)  *  dxx(i) 
20  CONTINUE 
C 

C      CALCULATE  HEAT  GENERATED  FROM  CHIP 
C 

do  30  m=l,njchip 
do  30  n=l, nkchip 

do  30  i=ibgn,iend 

do  30  j=jbgn (m) , jend (m) 
do  30  k=kbgn (n) ,kend (n) 

qchip  =  qchip  +  ale  *  dxx(i)  *  dyy ( j )  *  dzz(k) 
30   continue 

AA  =  1.  /  (bth  *  wth) 
RNUH  =  RNUH  *  AA 
RNUC  =  RNUC  *  AA 
QCHIP  =  QCHIP  *  AA 
QCRT  =  QCRT  *  AA 
QALL  =  QCHIP  +  RNUH  -  RNUC 

WRITE(6,500)  NT,XTIME, ITER, RESORM ( ITER) ,SORSUM, 
&  RNUC, RNUH, QCHIP, QALL, QUNS, QCRT 

500  FORMAT(lX,  'NT  =  ' , 19 , 5X, 'TIME=' , f 7 . 4 , / , 

&         IX,  'ITER=',I2,   5X, 'SOURCE=' ,F9.6,5X, 'SORSUM=' , F9 . 6 , / , 
&         IX,  'NUC='fF10.6,  5X, 'NUH=' ,F10.6,/, 
&         IX,  'QCHIP=' ,F10.6,5X, 'QALL=' ,F10.6,/, 
&         IX,  'QUNS=' ,F10.6, 5X, 'QCRT=' ,F10.6, /) 
END 

SUBROUTINE  SIP ( 1ST , JST , KST , ISP , JSP , KSP , PHI ) 

implicit  real*8  (a-h,o-z) 

common/ limits/ni , nipl , niml , nj , njpl , njml, nk, nkpl, nkml , 
&  nip2 ,njp2 ,nkp2 , iter ,nnmax 

common / to 1/ small , eps, sormax 

common /dims/h, wth, bth , hchip, wchip, bchip,bsub 

common/ coeff /ap(12, 16, 16)  ,aw(12,16,16)  ,ae(12,16,16)  , 
&  as (12, 16, 16) ,an(12, 16, 16) , ab ( 12 , 16 , 16) , 
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&  af (12,16,16) ,su(12, 16, 16) 

common/ abc/kmax , nl , kl 

common/count/nt,nmax, imax, itmax,krun, nprint 

COMMON/  MAT/A1(3072) ,A2(3072) , A3 (3072) , A4 ( 3072 ), A5 ( 3072 ) , 
&  A6 (3072) ,A7 (3072) ,B(3072) ,D (3072) 

DIMENSION  Xl(3072) ,X2(3072) ,XR1(3072) , 
&  S(12,16,16) ,PHI(12,16,16) 

C 

Nl   =  ISP  -  1ST  +  1 

Kl   =  (ISP  -  1ST  +  1)  *  (JSP  -  JST  +  1) 

KMAX  =  (ISP  -  1ST  +  1)  *  (JSP  -  JST  +  1)  *  (KSP  -  KST  +  1) 

EPS2  =  KMAX  *  EPS**2 
C 
C    DEFINING  A  SCRATCH  ARRAY  FOR  FURTHER  MODIFICATIONS 

C 

DO  9  K=KST,KSP 

DO  9  J=JST,JSP 

DO  9  I=IST,ISP 

S(I,J,K)  =  SU(I,J,K) 
9  CONTINUE 
C 

C    EXTRACT  MATRIX  A,  RHS  VECTOR  B  AND  VECTOR  XU-AN  INITIAL 
C  GUESS  FROM  GIVEN  DATA 

C 

DO  2  K=KST,KSP 

DO  2  J=JST,JSP 

DO  2  I=IST,ISP 

M   =  I  -  1ST  +  1  +  (J  -  JST)  *  Nl  +  (K  -  KST)  *  Kl 
A1(M)  =  AP(I,J,K) 
A2(M)  =  -AE(I,J,K) 
A3(M)  =  -AW(I,J,K) 
A4(M)  =  -AN(I,J,K) 
A5(M)  =  -AS(I,J,K) 
A6(M)  =  -AF(I,J,K) 
A7(M)  =  -AB(I,J,K) 
B(M)  =  S(I,J,K) 
X1(M)  =  PHI(I,J,K) 
2  CONTINUE 


c        THE  DIAGONAL  ELEMENT  OF  L  IN  THE  INCOMPLETE  LU 

c  DECOMPOSITION 

c 

D(l)  =  Al(l) 

DO  113  1=2, Nl 

D(I)  =  A1(I)  -  A2(I-1)  *  A3(I)  /  D(I-l) 

113  CONTINUE 

DO  114  I=N1+1,K1 

D(I)  =  A1(I)  -  A2(I-1)  *  A3(I)  /  D(I-l)  -  A4(I-N1) 
&  *  A5(I)  /  D(I-N1) 

114  CONTINUE 

DO  115  I=K1+1,KMAX 

D(I)  =  A1(I)  -  A2(I-1)  *  A3(I)  /  D(I-l) 
&  -  A4(I-N1)  *  A5(I)  /  D(I-N1) 

&  -  A6(I-K1)  *  A7(I)  /  D(I-K1) 

115  CONTINUE 

DO  116  1=1, KMAX 

D(I)  =  1.  /  D(I) 

116  CONTINUE 
ITR  =  0 
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c  THE  MAIN  ITERATION  LOOP  BEGINS 

c 

7   CONTINUE 

ITR  =  ITR  +  1 
CALL  RES  (X1,XR1) 
CALL  XL  (XR1,X2) 
DO  117  I=1,KMAX 

X2(I)  =  X2(I)  +  X1(I) 
117  CONTINUE 
V2  =  0. 
DO  55  I=1,KMAX 

V2  =  V2  +  XR1(I) **2 

55  CONTINUE 
X2MAG  =  0. 

DO  56  I=1,KMAX 

X2MAG  =  X2MAG  +  X2(I)**2 

56  continue 

if  (x2mag  .gt.  0.0001)  v2  =  v2  /  x2mag 

c 

c  RELABELING  THE  VECTORS  BEFORE  NEXT  ITERATION 

c 

do  50  i=l,kmax 
X1(I)  =  X2(I) 
50   continue 

IF  (V2  .GT.  EPS2  .AND.  ITR  .LT.  IMAX)  GO  TO  7 
c 

c     RECOVER  THE  SOLUTION  VECTOR  AND  RETURN  THE  AVERAGE 
c  PHI  AS  WELL  AS  THE  RMS  ERROR 

c 

IF  (MOD(NT,nprint)  .EQ.  0)  PRINT  *,  ITR,V2 
do  3  k=kst,ksp 
do  3  j=jst,jsp 
do  3  i=ist,isp 

m   =  i  -  ist  +  1  +  (j  -  jst)  *  nl  +  (k  -  kst)  *  kl 
phi(i, j,k)  =  xl(m) 
3  continue 
end 

SUBROUTINE  RES(X1,X2) 
implicit  real*8  (a-h,o-z) 
COMMON/ABC/KMAX , Nl , Kl 

common  / limits/ ni , nipl ,  niml,nj,njpl,njml,nk,  nkpl,  nkml, 
&  nip2 , njp2 , nkp2 , iter , nnmax 

common/  mat/al(3072) , a2 ( 3  072) , a3 ( 3072 )  , a4 ( 3072 )  , a5 (3072 ) , 
&  a6 (3  072) ,a7( 3  072) ,b(3072) ,d (3072) 

dimension  xl (3072) , x2 (3072) 
C 
c 

c      CALCULATE  Axl 
c 

X2(l)  =  al(l)  *  xl(l)  +  a2(l)  *  xl(2) 
&       +  a4(l)  *  xl(l+nl)  +  a6(l)  *  xl(l+kl) 
do  10  i=2,nl 

x2(i)  =  al(i)  *  xl(i)  +  a2(i)  *  xl(i+l) 
&         +  a4(i)  *  xl(i+nl)  +  a6(i)  *  xl(i+kl) 
&         +  a3(i)  *  xl(i-l) 
10   continue 

do  20  i=nl+l,kl 

x2(i)  =  al(i)  *  xl(i)  +  a2(i)  *  xl(i+l) 
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&         +  a4(i)  *  xl(i+nl)  +  a5(i)  *  xl(i-nl) 
&         +  a3(i)  *  xl(i-l)  +  a6(i)  *  xl(i+kl) 
20   continue 

DO  30  I=K1+1,KMAX-K1 

x2(i)  =  al(i)  *  xl(i)  +  a2(i)  *  xl(i+l) 
&         +  a4(i)  *  xl(i+nl)  +  a5(i)  *  xl(i-nl) 
&         +  a3(i)  *  xl(i-l)  +  a6(i)  *  xl(i+kl) 
&         +  a7(i)  *  xl(i-kl) 

30  continue 

DO  32  I=KMAX-K1+1,KMAX-N1 

x2(i)  =  al(i)  *  xl(i)  +  a2(i)  *  xl(i+l) 
&         +  a4(i)  *  xl(i+nl)  +  a5(i)  *  xl(i-nl) 
&  +  A3(I)  *  Xl(I-l)  +  A7(I)  *  Xl(I-Kl) 

3  2   CONTINUE 

DO  34  I=KMAX-N1+1,KMAX-1 

x2(i)  =  al(i)  *  xl(i)  +  a2(i)  *  xl(i+l) 
&  +  A5(I)  *  Xl(I-Nl)  +  A3(I)  *  Xl(I-l) 

&  +  A7(I)  *  Xl(I-Kl) 

3  4   CONTINUE 

X2(KMAX)  =  Al(KMAX)  *  Xl(KMAX)  +  A5(KMAX)  *  Xl(KMAX-Nl) 
&  +  A3(KMAX)  *  Xl(KMAX-l)  +  A7(KMAX)  *  Xl(KMAX-Kl) 

do  40  i=l,kmax 

X2(I)  =  B(I)  -  X2(I) 
40   continue 
end 
c* **************************************************************************** 

subroutine  xl(xl,x3) 
implicit  real*8  (a-h,o-z) 
common/abc/kmax, nl ,kl 

common/  mat/al(3072)  , a2 ( 3  07  2  )  , a3  (  3  07  2 )  , a4 ( 3  07  2 )  , a5 ( 3  072 )  , 
&  a6  (3072)  ,a7 (3072)  ,b(3072)  ,d (3072) 

dimension  xl(3072) ,x2(3072) ,x3(3072) 
X2(l)  =  xl(l)  *  d(l) 
do  31  i=2,nl 

X2(I)  =  D(I)  *  (X1(I)  -  A3(I)  *  X2(I-1)) 

31  continue 

do  32  i=nl+l,kl 

X2(I)  =  D(I)  *  (X1(I)  -  A3(I)  *  X2(I-1) 
&  -  A5(I)  *  X2(I-N1) ) 

32  continue 

do  33  i=kl+l,kmax 

X2(I)  =  D(I)  *  (X1(I)  -  A3(I)  *  X2(I-1) 
&  -  A5(I)  *  X2(I-N1)  -  A7(I)  *  X2(I-K1)) 

33  continue 
x3(kmax)  =  x2(kmax) 

do  67  i=kmax-l ,kmax-nl+l , -1 

X3(I)  =  X2(I)  -  D(I)  *  A2(I)  *  X3(I+1) 

67  continue 

do  68  i=kmax-nl ,kmax-kl+l , -1 

X3(I)  =  X2(I)  -  D(I)*(A2(I)*X3(I+1)  +  A4 ( I) *X3 ( I+Nl) ) 

68  continue 

do  69  i=kmax-kl , 1 , -1 

X3(I)  =  X2(I)  -  D(I)*(A2(I)*X3(I+1)  +  A4 ( I ) *X3 ( I+Nl ) 
&  +  A6(I)*X3(I+K1) ) 

69  continue 
end 

c****************************************************************************' 

SUBROUTINE  PROP 
C 
C      IN  THIS  SUBROUTINE  THE  INTERFACIAL  THERMAL  DIFFUSIVITIES  ARE 
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CALCULATED. 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
,  dxxs(40) ,dyys(40) ,dzzs(40) 

common/ parm/r a, pr, sorsum, dtime, xper , roll,dt_inv, xtime 
common/tol/small , eps , sormax 

common/dims/h,wth,bth, hchip, wchip, bchip, bsub 
common/ limits/ni, nipl, niml , nj , njpl , njml, nk, nkpl , nkml , 
t  nip2 ,njp2 , nkp2 , iter , nnmax 

common/ count /nt, nmax, imax, itmax, krun, nprint 

common/chip/ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) , nchp, 
,  isub, ichip, jchip, kchip 

common/dif f /ale, a Is, thot, tcool, tavg 
COMMON/RHOCP/RHS , RHC 

common/condu/alpha(12, 16,16) ,rho(12, 16, 16) ,visco(12, 16, 16) 
COMMON/POWER/QQQ,QCOND 
common/array/njehip, nkchip, ichoice 


if  (ichoice  . eq .  1)  then 

cond  =  0.065  -  7.895E-5  *  tavg 

rh  =  (1.825  -  0.00246  *  tavg)  *  1000. 

cp  =  4186.  *  (0.2411  +  3.7037E-4  *  tavg) 

alp  =  cond  /  (rh  *  cp) 

beta  =  0.00246  /  (1.825 

PRINT  *  t  ' ********    FLUID 
endif 
if  (ichoice  . eq .  2)  then 

cond  =  0.071 

rh  =  (2.002  -  0.00224  * 

cp  =  4186.  *  (0.2411  +  : 

alp  =  cond  /  (rh  *  cp) 

beta  =  0.00224  /  (2.002  -  0.00224  *  tavg) 

PRINT  *,'********  FLUID  IS  FC71  ************' 
endif 


-  0.00246  *  tavg) 

IS  FC75  ************' 


tavg)  *  1000. 
1.7037E-4  *  tavg) 


if  (ichoice  .ne.  0)  then 

ra  =  9.81  *  beta  *  qqq  *  h**2  *  1.0E-6 
t       /  (hchip  *  bchip  *  wchip  *  ale  *  cond  *  alp  *  alp) 

qcond  =  qqq  *  1.0E3 
i  /  (h  *  hchip  *  bchip  *  wchip  *  ale  *  cond) 

PRINT  *,'THE  POWER  DISSIPATED  PER  CHIP  IS', QQQ 

PRINT  *,'THE  BOUSSINESQ  NUMBER  IS',RA 
endif 


PRINT  THE  INPUT  PARAMETERS 


if  (ichoice 
print  * 
print  * 
print  * 

endif 


eq.  0)  then 

'A  CONSTANT  VISCOSITY  FLUID  IS  ASSUMED' 

WHOSE  PRANDTL  NUMBER  IS',pr 

AND  WHOSE  BOUSSINESQ  NUMBER  IS',ra 


PRINT  *,'THE  ASPECT  RATIOS  ARE  ' , BTH , WTH 
PRINT  *,'THE  TIME  STEP  IS', DTIME 
PRINT  *, 'NUMBER  OF  TIME  STEPS', NMAX 
if  (krun  .eq.  1)  then 

print  *, '*********THIS  IS  A  RESTART*********' 
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else 

print  *, '*********THIS  IS  NOT  A  RESTART*****' 
endif 
if  (nchp  .eq.  1)  then 

print  *,'*****cHIP  AND  SUBSTRATE  PRESENT*****' 
else 

print  *,'****NO  CHIP  OR  SUBSTRATE***' 
endif 
C 

C  SET  UP  THE  THERMAL  DIFFUSIVITY 

C  FOR  DIFFERENT  CONTROL  VOLUMES 

C 

C  FLUID  DIFFUSIVITY 

C 

DO  10  I=1,NIP1 

DO  10  J=1,NJP1 

DO  10  K=1,NKP1 

ALPHA (I, J, K)  =  1. 
RHO(I,J,K)  =  1. 
10   CONTINUE 

if  (nchp  .ne.  0)  then 
C 

C  SUBSTRATE  DIFFUSIVITY 

C 

DO  15  I=1,IBGN-1 
DO  15  J=1,NJP1 

DO  15  K=1,NKP1 

ALPHA (I, J, K)  =  ALS 
RHO(I,J;K)  =  RHS 
15       CONTINUE 
C 

C  CHIP  DIFFUSIVITY  FOR  CHIPS 

C 

DO  20  M=1,NJCHIP 

DO  20  N=1,NKCHIP 

do  20  i=ibgn,iend 

do  20  j=jbgn(m) , jend (m) 
do  20  k=kbgn (n) ,kend (n) 
ALPHA (I, J ,K)  =  ALC 
RHO(I,J,K)  =  RHC 
20      continue 
endif 
end 

subroutine  openf 
c 

c      THIS  SUBROUTINE  OPENS  ALL  FILES 
c 

c      UNIT  8:  THIS  IS  THE  INPUT  FILE  FOR  THE  FIELD  VARIABLES.  IT  IS 
c  USED  FOR  A  RESTART  JOB. 

c      UNIT  10:  THIS  FILE  STORES  THE  U,V  AND  W  VELOCITIES  FOR  A  GIVEN 
c  LOCATION. 

c      UNIT  11:  THIS  FILE  STORES  THE  FIELD  VARIABLES  AT  PRESCRIBED 
c  INTERVALS.  THIS  IS  A  CHECKPOINTING  PROCEDURE, 

c      UNIT  12:  THIS  FILE  STORES  THE  MEAN  FIELD  VARIABLES, 
c 

OPEN ( 8 , FILE= ' c2 1 . data ' , STATUS= ' UNKNOWN ' , FORM= ' UNFORMATTED ' ) 

OPEN ( 10 , FILE= ' c2 v . data ' , STATUS= ' UNKNOWN ' ) 

OPEN (11, FILE= ' c2 tmp . data ' , STATUS= ' UNKNOWN ' , FORM= ' UNFORMATTED ' ) 

OPEN ( 12, FILE='c2m. data' , STATUS=' UNKNOWN' , FORM= ' UNFORMATTED ' ) 
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0PEN(13 , FILE='c2grd.data' , STATUS=' UNKNOWN' ) 
end 
c******************************************************* 

subroutine  initio 
c 

c      IN  THIS  SUBROUTINE  ALL  VARIABLES  ARE  INITIALIZED  TO  DEFAULTS 
c      OR  THE  FIELD  VARIABLES  ARE  READ. 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common /parm/r a , pr , sor sum, d time, xper , roll ,dt_inv, xtime 

common /tol /small, eps , sormax 

common/dims/h, wth, bth, hchip,wchip, bchip, bsub 

common/ limits/ni , nipl ,niml,nj,njpl,njml,nk, nkpl, nkml, 
&  nip2 , njp2 , nkp2 , iter , nnmax 

common/ fv_init/ tod (12, 16, 16) ,uod(12,16, 16) , vod(12, 16, 16) , 
&  wod(12, 16, 16) ,pod(12, 16, 16) 

common/ fv_cur/t( 12, 16, 16) ,u(12,16,16) ,v(12,16,16) , 
&  w(12, 16, 16) ,p(12, 16, 16) 

common/fv_int/tpd(12, 16, 16) ,upd(12, 16, 16) , vpd(12, 16, 16) , 
&  wpd(12,16,16) ,ppd(12, 16,16) 


common /mscn/smp( 12,16,16 
&  du(12,16,16) 

&  dw(12,16,16) 

common /coef f /ap( 12,16,16 
&  as(12, 16, 16 

&  af(12,16,16 

common/mean/ t_mean ( 12  ,  16 
&  v_mean(12 , 16 

&  p  mean (12,16 


, resorm (93 ) , 

dv(12, 16, 16) , 

pp(12,16,16) 

,aw(12,16,16) ,ae(12,16,16) , 

,  an (12, 16, 16), ab (12, 16, 16), 

,su(12, 16, 16) 

16) ,u_mean( 12 , 16, 16) , 

16) ,w_mean(12 , 16, 16)  , 

16) 

common / count /nt, nmax, imax, it max, krun, nprint 

common/chip/ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) , kend (3) , nchp, 
&  isub, ichip, jchip, kchip 

common/scheme/quick, upwind 

common/dif f /ale, a Is, thot , tcool , tavg 

common / arr ay /n jchip, nkchip, i choice 

common/ if irst/nust, nvst , nwst , npst 

DATA  PI/3.14159265359/ 

SETTING  UP  THE  WEIGHTED  AVERAGE 

UPWIND  =  1.0  -  QUICK 

qck  =  100.  *  quick 

uwd  =  100.  *  upwind 

print  *,'THE  SCHEME  IS ', qck, ' PERCENTAGE  QUICK  AND' 

print  *; uwd, 'PERCENTAGE  UPWIND' 

SET  THE  I-DIRECTION  INDEX  FOR  CALU,CALV  AND  CALW 

then 


if  (nchp 

.eq.  0 

nust 

= 

3 

nvst 

= 

2 

nwst 

= 

2 

npst 

= 

2 

else 

nust 

= 

ibgn 

nvst 

= 

ibgn 

nwst 

= 

ibgn 

+  1 
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npst 
endif 


ibgn 


c 
c 
c 
c 


INITIALIZE  ALL  VARIABLES  TO  ZERO  OR  DEFAULTS 
BEFORE  START  OF  COMPUTATIONS 

DO  205  I=1,NIP1 

DO  205  J=1,NJP1 

DO  205  K=1,NKP1 

UOD(I,J,K)  =  0. 


U(I 
UPD 
VOD 
V(I 
VPD 
W(I 
WPD 
WOD 
POD 
P(I 
PPD 


=  0 


DU(I, J 
DV ( I , J 
DW ( I , J 
SU(I, J 
PP(I,J 
AP(I, J 
AW ( I , J 
AE(I, J 
AN ( I ,  J 
AS (I,  J 


J,K)  =  0 

I,J,K)  = 

I,J,K)  = 

J,K)  = 

I,J,K) 

J,K) 

I,J,K) 

I,J,K) 

I,J,K) 

J,K) 

I,J,K) 

K 

K 


=  0 
=  0 
0. 
=  0 


=  0 


=  0 


SMP(I, J,K) 


T_MEAN 
U_MEAN 
V_MEAN 
W_MEAN 
P  MEAN 


0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
=  0 


J,K) 
J,K) 
J,K) 
J,K) 
J,K) 


0. 
0. 
0. 
0. 
0. 


205  CONTINUE 


c 
C 

c 


RESTART  JOB  IF  KRUN  IS  SET  TO  ONE 

IF  (KRUN  . EQ.  1)  THEN 

READ (8)  TOD ,UOD, VOD, WOD, POD 
ELSE 

DO  220  J=1,NJP1 
DO  220  I=1,NIP1 

DO  220  K=1,NKP1 

TOD (I, J, K)  =  THOT  - 

*  (y(j)  - 

=  0. 

TOD(I, J,K) 
=  TOD (I, J, K) 


220 


tod(i, j,k) 
T(I,J,K)  = 
TPD(I, J,K) 
CONTINUE 


THOT  -  TCOOL) 
0.5  *  dyy(j)) 


c 
c 

c 


AMPLITUDE  OF  PERTURBATIONS 

Al  =  1.  /  bth 

YPER  =  ROLL  *  XPER  /  bth 

U  VELOCITY  PERTURBATIONS 
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do  230  k=2,nk 

DO  23  0  J=2,NJ 

DO  230  I=3,NI 

U(I,J,K)  =  -  xper  *  cos(p 

ii 

*  (y(j)  +  0.5  *  dyy(j)) ) 

&                   *  sin(roll  *  pi 

* 

al  *  x(i) ) 

UOD(I,J,K)  -  U(I,J,K) 

UPD(I,J,K)  =  U(I,J,K) 

230 

CONTINUE 

c 
c 

V  VELOCITY  PERTURBATIONS 

c 

do  235  k=2,nk 

DO  2  35  J=3,NJ 

DO  235  I=2,NI 

V(I,J,K)  =  yper  *  sin(pi 

* 

y(j)) 

&                      *  COS (ROLL  *  PI 

* 

Al  *  (x(i)  +  0.5  *  dxx(i) ) ) 

VOD(I,J,K)  =  V(I,J,K) 

VPD(I,J,K)  =  V(I,J,K) 

235    CONTINUE 
ENDIF 
C 

c  SET  THE  VELOCITIES  OUTSIDE  THE  COMPUTATIONAL 

c  DOMAIN  TO  ZERO 

c 

do  505  i=l,nipl 

DO  505  J=1,NJP1 

UOD(I,J,l)  =  0. 
UOD(I, J,NKP1)  =  0. 
VOD ( I , J , 1 )  =  0 . 
VOD(I, J,NKP1)  =  0. 
WOD(I,J,l)  =  0. 
WOD(I,J,2)  =  0. 
WOD(I, J,NKP1)  =  0. 
505  CONTINUE 

DO  510  I=1,NIP1 

DO  510  K=1,NKP1 

UOD(I, 1,K)  =  0. 
U0D(I,NJP1,K)  =  0. 
VOD(I,l,K)  =  0. 
VOD (I, 2, K)  =  0. 
V0D(I,NJP1,K)  =  0. 
WOD(I, 1,K)  =  0. 
W0D(I,NJP1,K)  =  0. 
510  CONTINUE 

DO  520  J=1,NJP1 

DO  520  K=1,NKP1 

UOD(l,J,K)  =  0. 
UOD(2,J,K)  =  0. 
U0D(NIP1, J,K)  =  0. 
VOD(l,J,K)  =  0. 
V0D(NIP1, J,K)  =  0. 
WOD(l,J,K)  =  0. 
W0D(NIP1, J,K)  =  0. 
520  CONTINUE 
c 

c  SET  THE  VELOCITIES  AND  PRESSURE  IN  THE  SUBSTRATE  TO  ZERO 

c 

if  (nchp  .ne.  0)  then 
do  530  i=l , ibgn 
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do  530  j=l,njpl 
do  530  k=l,nkpl 
uod(i, j , k)  =  0 . 
530     continue 
c 

do  540  i=l,ibgn-l 
do  540  j=l,njpl 
do  540  k=l,nkpl 
vod (i, j ,k)  =  0 . 
wod( i , j ,k)  =  0 . 
pod(i, j,k)  =  0. 
540      continue 
endif 
C 

C         INITIALISE  U,V,W,T,P 
C 

DO  210  K=1,NKP1 

DO  210  J=1,NJP1 

DO  210  I=1,NIP1 

T(I,J,K)  =  TOD(I,J,K) 
U(I,J,K)  =  UOD(I,J,K) 
V(I,J,K)  =  VOD(I,J,K) 
W(I,J,K)  =  WOD(I,J,K) 
P(I,J,K)  =  P0D(I,J,K) 
210  CONTINUE 
end 

subroutine  ploop 
c 

c      THIS  SUBROUTINE  INCLUDES  THE  PRESSURE  LOOP.  IT  INCORPORATES  THE 
c      SIMPLEX  ALGORITHM  AND  THE  INFAMOUS  ERROR  CONTROL  ROUTINE  DUE  TO 
c      THE  VENERABLE  VINCENT  LIU  OF  THE  ARGONNE  NATIONAL  LABS  WHO  DEVISED 
c      IT  AS  A  HOBBY  WHILE  AT  THE  UNIV  OF  NOTRE  DAME. 


c 


implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40) , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/ limits/ni  ,nipl,niml,nj,njpl,njml,nk,  nkpl,nkml, 
&  nip2 , njp2 ,nkp2 , iter , nnmax 

common/ pa rm/ra, pr , sorsum, dtime, xper , roll,dt_inv,xtime 

common/ tol/ small , eps, sormax 

common/dims/h,wth,bth,hchip, wchip,bchip, bsub 

common/ fv_init/tod( 12, 16, 16) ,uod(12, 16, 16) , vod (12, 16, 16) 
&  wod(12, 16, 16) ,pod(12, 16, 16) 

common/ fv_cur/t( 12, 16,16),u(12,16,16),v(12,16,16), 
&  w(12, 16, 16) ,p(12, 16, 16) 

common/ fv_int/tpd( 12, 16, 16) ,upd(12, 16, 16) , vpd(12, 16, 16) , 
&  wpd(12, 16, 16) ,ppd(12, 16, 16) 

common /mscn/smp( 12 , 16, 16) , resorm(93 ) , 
&  du(12, 16, 16) ,dv(12, 16, 16) , 

&  dw(12, 16, 16) ,pp(12, 16, 16) 

common/ coef f /ap(12 ,16, 16) , aw (12, 16, 16) ,ae(12,16,16) , 
&  as(12,16,16),an(12,16,16),ab(12,16,16), 

&  af (12, 16, 16) ,su(12, 16, 16) 

COMMON/COEF2/AWW(12, 16, 16) ,AEE(12, 16, 16) , ASS (12, 16, 16) , 
&  ANN (12, 16,16) , ABB (12, 16, 16) ,AFF(12, 16, 16) 

common / mean /t_mean (12,16,16) , u_mean(12 , 16 , 16) , 
&  v_mean( 12 , 16 , 16) , w_mean (12,16,16) , 

&  p_mean ( 12 , 16, 16) 

common/ count /nt , nmax, imax, itmax,krun, nprint 
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common/chip/ ibgn,iend, jbgn(3) , jend(3) ,kbgn(3) , kend(3) ,nchp, 
&  isub, ichip, jchip,kchip 

common/array/njchip,nkchip, ichoice 
common/diff /ale, als, thot , tcool,tavg 

c 

DIMENSION  UU(5000)  ,W(5000)  ,WW(5000) 
c        CALL  CPUTIME  (BEGIN, IIR) 
XTIME  =  0. 
NT  =  0 
c 

c     THE  MAIN  TIME  MARCHING  LOOP  BEGINS  HERE 
c 

300  continue 
c 

NT  =  NT  +  1 
c 

c     UPDATE  THE  MEAN  VELOCITIES , TEMPERATURE  AND 
c  PRESSURE 

c 

CI  =  1.  /  FLOAT (NT) 
do  310  k=l,nkpl 

DO  310  J=1,NJP1 

DO  310  I=1,NIP1 

U_MEAN(I,J,K)  =  (1.  -  CI)  *  U_MEAN(I,J,K) 
&  +  CI  *  U(I,J,K) 

V_MEAN(I,J,K)  =  (1.  -  CI)  *  V_MEAN(I,J,K) 
&  +  cl  *  v(i, j,k) 

W_MEAN(I, J,K)  =  (1.  -  Cl)  *  W_MEAN(I,J,K) 
&  +  cl  *  w(i, j,k) 

T_MEAN(I,J,K)  =  (1.  -  Cl)  *  T_MEAN(I,J,K) 
&  +  cl  *  t(i, j ,k) 

310  continue 
c 

c       STORE  THE  FIELD  VARIABLES  EVERY  1000  TIME  STEPS 
c  AS  A  KIND  OF  CHECKPOINTING 

c 

IF  (MOD(NT,nprint)  .EQ.  0)  THEN 
WRITE(ll)  T,U,V,W,P 
REWIND  11 
ENDIF 
c 

C      STOP  COMPUTATIONS  AT  THE  MAXIMUM  PRESCRIBED 
c  TIME  STEP 

c 

IF(NT  .GT.  NMAX)  THEN 
DO  10  1=1, NMAX 

WRITE(10,*)  UU(I) ,VV(I) ,WW(I) 
10     CONTINUE 

WRITE ( 12 )  T_MEAN , U_MEAN , V_MEAN , W_MEAN , P_MEAN 
c  CALL  CPUTIME  (END, IIR) 

TT  =  (END  -  BEGIN)  *  l.E-06 
PRINT  *,'THE  CPUTIME  USED  WAS',TT 
STOP 
ENDIF 
C 

XTIME  =  XTIME  +  DTIME 
C 

C         START  CALCULATIONS 
C 

ITER  =  0 
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JTERM  =  0 
JJTERM  =  0 

DEFINE  THE  UPDATED  TPD(I,J,K),  UPD(I,J,K)  AND  VPD(I,J,K) 

FOR   SU(I,J,K) 

DO  48  K=1,NKP1 

DO  48  J=1,NJP1 

DO  48  I=1,NIP1 

TPD(I,J,K)  =  T(I,J,K) 
UPD(I,J,K)  =  U(I,J,K) 
VPD(I,J,K)  =  V(I,J,K) 
WPD(I,J,K)  =  W(I,J,K) 
48  CONTINUE 

29  CONTINUE 

JTERM  =  JTERM  +  1 

CALL  CALT 

DO  2220  I=1,NIP1 

DO  2220  K=1,NKP1 
DO  2  2  20  J=2,NJ 

IF  (T(I,J,K)  .LT.  TCOOL)  T(I,J,K)  =  TCOOL 
CONTINUE 

PRESSURE  CORRECTION  LOOP 

301  CONTINUE 

ITER  =  ITER  +  1 

CALL  CALU 

CALL  CALV 

CALL  CALW 

CALL  CALP 

if (resorm(iter)  .le.  sormax)  go  to  49 
if(iter  .eq.  1)  go  to  302 

if (resorm(iter)  .le.  resorm( iter-1) )  go  to  302 
go  to  304 

302  if(jterm  .ge.  2)  go  to  37 
source=resorm ( iter) 

go  to  39 

37  if (resorm(iter)  .le.  source)  go  to  38 
go  to  304 

38  source=resorm( iter) 

39  continue 

do  23  k=l,nkpl 

do  23  j=l,njpl 
do  23  i=l,nipl 

tpd(i,j,k)  =  t(i,j,k) 
upd(i, j,k)  =  u(i, j,k) 
vpd(i,j,k)  =  v(i,j,k) 
wpd(i, j ,k)  =  w(i, j,k) 
ppd(i, j,k)  =  p(i, j,k) 
23  continue 
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j jterm=0 

if(iter  .eq.  itmax)  go  to  49 
if(jterm  .eq.  2)  go  to  35 
if(iter  .eq.  4)  go  to  29 
35  continue 

if(jterm  .eq.  3)  go  to  58 
if(iter  .eq.  7)  go  to  29 
58  continue 
j jterm=0 
go  to  301 
304  continue 

j jterm=j jterm+1 
if(jterm  . eq .  1)  go  to  41 

if(jterm  . eq .  2  .and.  jjterm  . eq .  1  .and.  iter  .ne.  5)  go  to  41 
go  to  82 
41  continue 

do  40  k=l,nkpl 
do  40  j=l,njpl 

do  40  i=l,nipl 

u(i,j,k)  =  upd(i,j,k) 
v(i,j,k)  =  vpd(i,j,k) 
w(i, j,k)  =  wpd(i, j,k) 
P(i, j,k)  =  ppd(i,j,k) 
40  continue 

if(iter  .eq.  itmax)  go  to  49 
go  to  29 
82  continue 

do  43  k=l,nkpl 
do  43  j=l,njpl 
do  43  i=l,nipl 

t(i,j,k)  =  tpd(i,j,k) 
u(i, j,k)  =  upd(i, j ,k) 
v(i, j,k)  =  vpd(i, j ,k) 
w(i, j,k)  =  wpd(i, j,k) 
P(i, j,k)  =  PPd(i, j ,k) 
43  continue 

if(iter  .eq.  itmax)  go  to  49 

if((jterm  .eq.  3  .and.  iter  .ne.  8)  .or.  jjterm  .eq.  2)  go  to  49 
go  to  301 
49  continue 

II  =  NIP1  *  2  /  3 
JJ  =  NJP1  *  2  /  3 
KK  =  NKP1  *  2  /  3 
UU(NT)  =  U(II,JJ,KK) 
W(NT)  =  V(II/JJ,KK) 
WW(NT)  =  W(II,JJ,KK) 

PRINT  ENERGY  AND  MASS  BALANCE  STATISTICS  AT  REGULAR 

INTERVALS 

IF  (MOD(NT,nprint)  .EQ.  0)  THEN 
CALL  NU 
call  chiptemp 
call  tstep 
ENDIF 
c 

c         UPDATE  VARIABLES  FOR  THE  NEXT  TIME  STEP 
c 

DO  305  K=1,NKP1 
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DO  305  J=1,NJP1 

DO  305  I=1,NIP1 

TOD(I,J,K)  =  T(I,J,K) 
UOD(I,J,K)  =  U(I,J,K) 
VOD(I,J,K)  =  V(I,J,K) 
WOD(I,J,K)  =  W(I,J,K) 
305  CONTINUE 
GO  TO  300 
end 

subroutine  tstep 

implicit  real*8  (a-h,o-z) 
c 

c      THE  MAXIMUM  TIME  STEP  FOR  AN  EXPLICIT  MARCH  IS  CALCULATED, 
c      THIS  IS  USEFUL  FOR  DECIDING  THE  TIME  STEP  THAT  NEEDS 
c      TO  BE  SPECIFIED, 
c      ALSO,  THE  VOLUME  AVERAGED  PRANDTL  NUMBER  IS  CALCULATED. 


c 


common/ fv_cur/t( 12, 16,16) ,u ( 12 , 16 , 16) , v( 12 , 16 , 16) , 
&  w(12, 16, 16) ,p(12, 16, 16) 

common/condu/alpha(12, 16, 16) ,rho(12, 16, 16) , visco(12, 16, 16) 

common /dims/h, wth, bth,hchip, wchip, bchip, bsub 

common/ng/dxx(40) ,dyy(40) ,dzz(40),x(40),y(40),z(40), 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/ limits/ ni, nipl,niml,nj,njpl,njml,nk, nkpl, nkml, 
&  nip2 , njp2 , nkp2 , iter , nnmax 

common / pa rm/ r a, pr , sor sum, d time, xper , roll,dt_inv, xtime 

common/ if irst/nust , nvst , nwst , npst 


tmin  = 

-    0. 

do 

10 

i=2 , ni 

do 

10 

j=2,nj 

do 

10  k= 

=  2, 

,  nk 

tmpl 

= 

0.5 

* 

& 

+ 

& 

4- 

((u(i,j,k)  +  u(i+l,j,k))  /  dxx(i) 
(v(i, j,k)  +  v(i, j+i,k) )  /  dyy(j) 
(w(i,j,k)  +  w(i,j,k+l))  /  dzz(k)) 
tmp2    1.  /  dxx(i)**2  +  1.  /  dyy(j)**2 
&  +  1.  /  dzz(k)**2 

tmp  =  max  (tmpl,tmp2) 
tmp  =  tmpl 

if  (tmp  .gt.  tmin)  tmin  =  tmp 
10   continue 

tmin  =1.  /  tmin 

umin  =  0. 
do  20  i=nust,ni 
do  20  j=2,nj 

do  20  k=2,nk 

tmpl  =  0.5  *  (2.  *  u(i,j,k)  /  dxxs(i) 
&  +  (v(i,j,k)  +  v(i,j+l,k))  /  dyy(j) 

&  +  (w(i,j,k)  +  w(i,j,k+l))  /  dzz(k)) 

tmp2  =  visco(i,j,k)  *  (1.  /  dxxs(i)**2 
&  +  1.  /  dyy(j)**2  +  1.  /  dzz(k)**2) 

tmp  =  max  (tmpl,tmp2) 
tmp  =  tmpl 

if  (tmp  .gt.  umin)  umin  =  tmp 
20   continue 

umin  =  1.  /  umin 

vmin  =  0. 
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do  30  i=nvst,ni 
do  30  j=3,nj 

do  3  0  k=2,nk 

tmpl  =  0.5  *  (2.  *  v(i,j,k)  /  dyys(j) 
&  +  (u(i,j,k)  +  u(i+l,j,k))  /  dxx(i) 

&  +  (w(i,j,k)  +  w(i,j,k+l))  /  dzz(k)) 

tmp2  =  visco(i,j,k)  *  (1.  /  dxx(i)**2 
&  +  1.  /  dyys(j)**2  +  1.  /  dzz(k)**2) 

c  tmp  =  max  (tmpl,tmp2) 

tmp  =  tmpl 

if  (tmp  .gt.  vmin)  vmin  =  tmp 
30   continue 

vmin  =  1.  /  vmin 
c 

wmin  =0. 
do  40  i=nwst,ni 
do  40  j=2,nj 

do  40  k=3,nk 

tmpl  =  0.5  *  (2.  *  w(i,j,k)  /  dzzs(k) 
&  +  (v(i,j,k)  +  v(i,j+l,k))  /  dyy(j) 

&  +  (u(i,j,k)  +  u(i+l,j,k))  /  dxx(i)) 

tmp2  =  visco(i,j,k)  *  (1.  /  dxx(i)**2 
&  +  I-  /  dyy(j)**2  +  1.  /  dzzs(k)**2) 

c  tmp  =  max  (tmpl,tmp2) 

tmp  =  tmpl 

if  (tmp  .gt.  wmin)  wmin  =  tmp 
40   continue 

wmin  =  1.  /  wmin 
tt  =  min  (tmin,umin, vmin, wmin) 
print  *,'THE  EXPLICIT  TIME  STEP  IS',tt 
c 
c      THE  VOLUME  AVERAGED  PRANDTL  NUMBER  IS  CALCULATED. 

c 

vavgpr  =  0 . 
do  50  i=nust,ni 
do  50  3=2, nj 

do  50  k=2,nk 

vavgpr  =  vavgpr  +  dxx(i)  *  dyy(j)  *  dzz(k) 
&  *  visco( i , j ,k) 

50   continue 

vavgpr  =  vavgpr  /  (wth  *  bth) 

print  *,'THE  VOLUME  AVERAGED  PRANDTL  NUMBER  IS', vavgpr 

end 

c* ********************************************************************** 

subroutine  grid 
c 

C      THIS  SUBROUTINE  GENERATES  THE  GRID.  THE  ROBERTS  TRANSFORMATION 
c      IS  USED  TO  RESOLVE  BOUNDARY  LAYER  NEAR  THE  WALL. 


c 


implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40),x(40),y(40),z(40), 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common /blayer/xbr ,ybr , zbr 

common /unf rm/iunf rm 

common/dims /h, wth, bth, hchip, wchip, bchip, bsub 

common/chip/ibgn, iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) ,nchp, 
&  isub, ichip, jchip,kchip 

common/ limits/ni, nipl,niml,nj,njpl,njml,nk, nkpl, nkml, 
&  nip2 , njp2 ,nkp2 

common/ arr ay /nj chip, nkchip 
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common/spaced/ychip(3)  ,zchip(3) 

dimension  kgap(4) , jgap(4) ,ygap(4) ,zgap(4) , 
&  dygap(4) ,dzgap(4) 

print  *,'THE  HEIGHT  OF  THE  ENCLOSURE  IS',h,'mm' 
print  *,'THE  WIDTH  OF  THE  ENCLOSURE  IS ' , wth , ' mm" 
print  *,'THE  LENGTH  OF  THE  ENCLOSURE  IS', bth, 'mm' 

wth  =  wth  /  h 
bth  =  bth  /  h 
xbr  =  xbr  /  h 
ybr  =  ybr  /  h 
zbr  =  zbr  /  h 

print  *,'THE  CHIP  DIMENSIONS  ARE  THE  FOLLOWING:' 

print  *,'CHIP  HEIGHT  ( Y-DIRECTION) ' , hchip, 'mm' 

print  *,'CHIP  WIDTH  (Z-DIRECTION) ', wchip, 'mm' 

print  *,'CHIP  LENGTH (X-DIRECTION) ', bchip, 'mm' 

NONDIMENSIONALIZE  THE  LENGTH  PARAMETERS  WHICH  ARE  GIVEN  IN  MM 

hchip  =  hchip  /  h 
wchip  =  wchip  /  h 
bchip  =  bchip  /  h 
bsub  =  bsub  /  h 

if  (iunfrm  .eq.  0)  then 
if  (nchp  .ne.  0)  then 

do  10  i=l,njchip 

ychip(i)  =  ychip(i)  /  h 
continue 

do  20  i=l,nkchip 

zchip(i)  =  zchip(i)  /  h 
20   continue 

dx  =  bchip  /  float(ichip) 
dy  =  hchip  /  float(jchip) 
dz  =  wchip  /  float(kchip) 

zgap(l)  =  zchip(l)  -  0.5  *  wchip 

zgap(nkchip+l)  =  wth  -  zchip(nkchip)  -  0.5  *  wchip 

ygap(l)  =  ychip(l)  -  0.5  *  hchip 

ygap(njchip+l)  =  1.0  -  ychip(njchip)  -  0.5  *  hchip 

kgap(l)  =  int(f loat(kchip)  *  zgap(l)  /  wchip)  -  4 
dzgap(l)  =  zgap(l)  /  f loat (kgap( 1) ) 
jgap(l)  =  int(f loat( jchip)  *  ygap(l)  /  hchip)  +  3 
dygap(l)  =  ygap(l)  /  f loat ( jgap(l) ) 

do  25  i=2,nkchip 

zgap(i)  -  zchip(i)  -  zchip(i-l)  -  wchip 
kgap(i)  =  int (float (kchip)  *  zgap(i)  /  wchip)  -  2 
dzgap(i)  =  zgap(i)  /  float (kgap( i) ) 
25   continue 

do  30  i=2,njchip 

ygap(i)    =   ychip(i)    -   ychip(i-l)    -   hchip 

jgap(i)    =    int(f loat(jchip)    *    ygap(i)    /    hchip)    +    1 
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c 


c 


c 


dygap(i)  =  ygap(i)  /  f loat( jgap(i) ) 
30   continue 

ktotal  =  0 

jtotal  =  0 

do  35  i=l,nkchip 

ktotal  =  ktotal  +  kgap(i) 
35   continue 

do  40  i=l,njchip 

jtotal  =  jtotal  +  jgap(i) 
40   continue 

ktotal  =  ktotal  +  nkchip  *  kchip 

jtotal  =  jtotal  +  njchip  *  jchip 

kgap(nkchip+l)  =  nkml  -  ktotal 

dzgap(nkchip+l)  =  zgap(nkchip+l)  /  float (kgap(nkchip+l) ) 

jgap(njchip+l)  =  njml  -  jtotal 

dygap(njchip+l)  =  ygap(njchip+l)  /  float (jgap (n^chip+l) ) 

if  ((nkml  -  ktotal)  .It.  4)  then 

print  *, 'INSUFFICIENT  #  OF  POINTS  IN  THE  Z-DIRECTION' 

stop 
endif 

if  ((njml  -  jtotal)  .It.  6)  then 

print  *, 'INSUFFICIENT  #  OF  POINTS  IN  THE  Y-DIRECTION' 

stop 
endif 

call  gridl  (ygap(l) ,ybr, jgap(l) ,ykl,yhl) 

call  gridl  (ygap(njchip+l) , ybr , jgap(njchip+l) , yk2 , yh2) 

call  gridl  (zgap(l) , zbr , kgap ( 1) , zkl , zhl) 

call  gridl  (zgap(nkchip+l) , zbr , kgap(nkchip+l) , zk2 , zh2) 

GRID  GENERATION  IN  THE  Y  AND  Z-DIRECTION  (TOP  AND  BOTTOM) 

do  45  j=2, jgap(l)+2 

y(j)  =  tanh  (ykl  *  yhl  *  dygap(l)  *  float(j  -  jgap(l)  -  2)) 
&         /  ykl  +  ygap(l) 
45   continue 

y(D  =  "  y(3) 

do  50  k=2,kgap(l)+2 

z(k)  =  tanh  (zkl  *  zhl  *  dzgap(l)  *  float(k  -  kgap(l)  -  2)) 
&        /  zkl  +  zgap(l) 
50   continue 

z(l)  =  -  z(3) 

jj  =  2  +  jtotal 

yy  =  ychip(njchip)  +  0.5  *  hchip 

do  55  j=l, jgap(njchip+l) +1 

y(j+jj)  =  tanh  (yk2  *  yh2  *  dygap(njchip+l) 
&  *  float(j) )  /  yk2  +  yy 

55   continue 

y(njp2)  =  1.0  +  y(njpl)  -  y(nj) 

kk  =  2  +  ktotal 

zz  =  zchip(nkchip)  +  0.5  *  wchip 

do  60  k=l,kgap(nkchip+l) +1 

z(k+kk)  =  tanh  (zk2  *  zh2  *  dzgap (nkchip+1) 
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c 


c 


&  *  float(k) )  /  zk2  +  zz 

60   continue 

z(nkp2)  ■  wth  +  z(nkpl)  -  z(nk) 

yy  =  -  hchip 
jj  =  2  -  jchip 
do  65  i=l,njchip 

jj  *  jj  +  jchip  +  jgap(i) 
yy  =  yy  +  hchip  +  ygap(i) 
do  65  j=l, jchip 

y(jj+j)  =  dy  *  float(j)  +  yy 
65   continue 

zz  =  -  wchip 
kk  =  2  -  kchip 
do  70  i=l,nkchip 

kk  =  kk  +  kchip  +  kgap(i) 
zz  =  zz  +  wchip  +  zgap(i) 
do  70  k=l, kchip 

z(kk+k)  =  dz  *  float(k)  +  zz 
70   continue 

yy  =  0. 

jj  -  2 

do  75  i=2,njchip 

jj  =  jj  +  jchip  +  jgap(i-l) 
yy  =  yy  +  hchip  +  ygap(i-l) 
do  75  j=l,jgap(i) 

y(jj+j)  =  dygap(i)  *  float(j)  +  yy 
75   continue 

zz  =  0. 

kk  =  2 

do  80  i=2,nkchip 

kk  =  kk  +  kchip  +  kgap(i-l) 
zz  =  zz  +  zgap(i-l)  +  wchip 
do  80  k=l,kgap(i) 

z(kk+k)  =  dzgap(i)  *  float(k)  +  zz 

80   continue 
c 
c      GRID  GENERATION  IN  THE  X-DIRECTION 

c 

ii  =  (niml  -  ichip  -  isub)  /  2 

ii  =  2  *  ii 

inc  =  niml  -  ichip  -isub  -  ii 

isub  =  isub  +  inc 

dsub  =  bsub  /  float (isub) 

do  160  i=l,isub+2 

x(i)  =  dsub  *  float(i  -  2) 
160  continue 
c 

dx  =  bchip  /  float(ichip) 

xx  =  bsub 

ii  =  isub  +  2 

do  170  i=l, ichip 

x(i+ii)  =  dx  *  float(i)  +  xx 
170  continue 
c 

ii  =  isub  +  ichip  +  2 

bb  =  0.5  *  (bth  -  bsub  -  bchip) 


c 
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xx  =  bsub  +  bchip  +  bb 

iil  =  (niml  -  ichip  -  isub)  /  2 

call  gridl  (bb, xbr , iil , xk2 , xh2) 

ii2  =  2  *  iil 

dr  =  bb  /  float(iil) 

do  180  i=l,ii2 

x(i+ii)  =  tanh  (xh2  *  xk2  *  dr  *  float  (i  -  iil))  /  xk2 
&  +  xx 

180  continue 

x(nipl+l)  =  bth  +  x(nipl)  -  x(ni) 
c 

ibgn  =  isub  +  2 
iend  =  isub  +  ichip  +  1 
c 

jbgn(l)  =  jgap(l)  +  2 

jend(l)  =  jbgn(l)  +  jchip  -  1 

do  200  i=2,njchip 

jbgn(i)  =  jbgn(i-l)  +  jchip  +  jgap(i) 
jend(i)  =  jbgn(i)  +  jchip  -  1 
200  continue 
c 

kbgn(l)  =  kgap(l)  +  2 
kend(l)  =  kbgn(l)  +  kchip  -  1 
do  220  i=2,nkchip 

kbgn(i)  =  kbgn(i-l)  +  kchip  +  kgap(i) 
kend(i)  =  kbgn(i)  +  kchip  -  1 
220  continue 
c 

else 

print  *,'THE  HEIGHT  OF  THE  ENCLOSURE  IS',h,'mm' 
print  *,'THE  WIDTH  OF  THE  ENCLOSURE  IS ' , wth , 'mm' 
print  *,'THE  LENGTH  OF  THE  ENCLOSURE  IS', bth, 'mm' 

c 

c      NONDIMENSIONALIZE  THE  LENGTH  PARAMETERS  WHICH  ARE  GIVEN  IN  MM 

c 

c 

c      GRID  GENERATION  IN  THE  Z-DIRECTION 

c 

zl  =  wth  *  0.5 

kzl  =  (nkpl  -  2)  /  2 

drl  =  zl  /  float(kzl) 

call  gridl  ( zl , zbr , kzl , zkl , zhl) 
c 

z(l)  =  -  zbr 

ii  =  2  *  kzl  +  l 

do  181  i=l,ii 

z(i+l)  =  tanh  (zhl  *  zkl  *  drl  *  float  (i  -  kzl  -  1) )  /  zkl 
&  +  zl 

181  continue 

z(nkpl+l)  =  wth  +  zbr 
c 

c      GRID  GENERATION  IN  THE  Y-DIRECTION 
c 

yl  =  0.5 

jyl  =  (njpl  -  2)  /  2 

drl  =  yl  /  float(jyl) 

call  gridl  (yl , ybr , jyl , ykl , yhl) 
c 

y(l)  =  -  ybr 

ii  =  2  *  jyl  +  1 
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c 
c 
c 


do  182  i=l,ii 

y(i+l)  =  tanh  (yhl  *  ykl  *  drl  *  float  (i 

&  +  yl 

182  continue 

y(njpl+l)  ■  1.0  +  ybr 

GRID  GENERATION  IN  THE  X-DIRECTION 


xl  =  0.5  *  bth 

ixl  =  (nipl  -  2)  /  2 

drl  =  xl  /  float(ixl) 

call  gridl  (xl, xbr , ixl , xkl, xhl) 

x(l)  -  -  xbr 

ii  =  2  *  ixl  +  1 

do  183  i=l,ii 

x(i+l)  =  tanh  (xhl  *  xkl  *  drl  *  float  (i  -  ixl 
&  +  xl 

183  continue 

x(nipl+l)  =  bth  +  xbr 

endif 

else 

print  *,'THE  GRID  IS  UNIFORM' 

dx  =  bth  /  float(niml) 
dy  =  1.0  /  float (njml) 
dz  =  wth  /  float(nkml) 

UNIFORM  GRID 

do  185  i=l,nip2 

x(i)  =  dx  *  float(i  -  2) 
185  continue 


jyl  -  1))  /  ykl 


-  D)  /  xkl 


do  186  j=l,njp2 

y(j)  =  dy  *  float(j  -  2) 
186  continue 


do  187  k=l,nkp2 

z(k)  =  dz  *  float(k  -  2) 
187  continue 

IF  (NCHP  .NE.  0)  THEN 
IBGN  =  4 
IEND  =  5 


JBGN(1 
JEND(1 
JBGN(2 
JEND(2 
JBGN(3 
JEND(3 
KBGN(1 
KEND(1 
KBGN(2 
KEND(2 
KBGN ( 3 
KEND(3 

ENDIF 

endif 


=  4 
=  5 
=  8 
=  9 
=  12 
=  13 
=  4 
=  5 
=  8 
=  9 
=  12 
=  13 


CALCULATE  THE  DIMENSIONS  OF  THE  CONTROL  VOLUMES 
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c 

do    188    i=l,nipl 

dxx(i)    =   x(i+l)    -   x(i) 

188  continue 

do    189    i=l,njpl 

dyy(i)    =  y(i+i)    -   y(i) 

189  continue 

do  191  i=l,nkpl 

dzz(i)  =  z(i+l)  -  z(i) 
191  continue 

dxx(nip2)  -  dxx(nipl) 
dyy(njp2)  =  dyy(njpl) 
dzz(nkp2)  =  dzz(nkpl) 


CALCULATE  THE  DIMENSIONS  OF  THE  STAGGERED  CONTROL  VOLUMES 


do  245  i=2,nipl 

dxxs(i)  =  0.5  *  (dxx(i-l)  +  dxx(i)) 
245  continue 

dxxs(l)  =  dxxs(2) 
do  250  i=2,njpl 

dyys(i)  =  0.5  *  (dyy(i-l)  +  dyy(i)) 
250  continue 

dyys(l)  =  dyys(2) 
do  255  i=2,nkpl 

dzzs(i)  =  0.5  *  (dzz(i-l)  +  dzz(i)) 
255  continue 

dzzs(l)  =  dzzs(2) 
dxxs(nip2)  =  dxxs(nipl) 
dyys(njp2)  =  dyys(njpl) 
dzzs(nkp2)  =  dzzs(nkpl) 
c 

c      CALCULATE  THE  AREAS  AND  VOLUMES  TO 
c      CHECK  ACCURACY  OF  GRID 


c 


svol  =  0. 

do  260  i=2,nipl-l 

do  260  j=2,njpl-l 
do  260  k=2,nkpl-l 

svol  =  svol  +  dxx(i)  *  dyy(j)  *  dzz(k) 

260  continue 

svol  =  svol  *  h**3 

print  *,'THE  TOTAL  VOLUME  IS ', svol ,' cubic  mm' 
svol  =  0. 
do  261  i=2,nipl-l 
do  261  j=2,njpl-l 

svol  =  svol  +  dxx(i)  *  dyy(j) 

261  continue 

svol  =  svol  *  h**2 

print  *,'THE  TOTAL  XY  AREA  IS ', svol ,' square  mm' 

svol  =  0. 

do  262  j=2,njpl-l 

do  262  k=2,nkpl-l 

svol  =  svol  +  dyy(j)  *  dzz(k) 

262  continue 

svol  =  svol  *  h**2 

print  *,'THE  TOTAL  YZ  AREA  IS ', svol ,' square  mm' 
svol  =  0. 
do  263  k=2,nkpl-l 
do  263  i=2,nipl-l 
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svol  =  svol  +  dzz(k)  *  dxx(i) 
263  continue 

svol  =  svol  *  h**2 

print  *,'THE  TOTAL  ZX  AREA  IS ', svol ,' square  mm' 

INDIC  =  MAXO  (NIP2,NJP2,NKP2) 

DO  270  1=1, INDIC 

write (6, 1000)  i,dxx(i) ,x(i) ,dyy(i) ,y(i),dzz(i),z(i) 

270  continue 
1000  format(lx,i5,2x,6(f9.5, lx) ) 

PRINT  * , IBGN , IEND , JBGN , JEND , KBGN , KEND 

write(13,2000)  x,y,z,dxx,dyy,dzz 
2000  format(4(fl7.7) ) 

end 

subroutine  gridl  (bb,dx, ngrid, xl ,hk) 
implicit  real*8  (a-h,o-z) 

C      THIS  SUBROUTINE  CALCULATES  THE  TWO  GRID  PARAMETERS 

c 

data  epps,nmax/1.0e-8, 1000/ 

dr  =  bb  /  float(ngrid) 

c  =  1.  /  (float(nmax)  *  bb) 

yy  =  1.  -  dr  /  bb 

yl  =  yy  -  1. 
c 
C      NARROW  DOWN  THE  ZERO  USING  THE  BISECTION  METHOD 

c 

i  =  0 
ii  =  0 
8   continue 
i  =  i  +  1 
xl  =  c  *  float(i) 

fl  =  (1.  +  xl  *  (bb  -  dx) )  *  (1.  -  xl  *  bb)**yy 
&    -  (l.  -  xl  *  (bb  -  dx) )  *  (1.  +  xl  *  bb)**yy 
x2  =  c  *  float(i+l) 

f2  =  (1.  +  x2  *  (bb  -  dx) )  *  (1.  -  x2  *  bb)**yy 
&    -  (l.  -  x2  *  (bb  -  dx))  *  (1.  +  x2  *  bb)**yy 
if  (fl  *  f2  .It.  0.  .or.  i  .eq.  (nmax  -  1))  then 

go  to  9 
endif 
go  to  8 
9    continue 

x3   =  0.5  *  (xl  +  x2) 
ii  =  ii  +  1 

f3  =  (1.  +  x3  *  (bb  -  dx) )  *  (1.  -  x3  *  bb) **yy 
&    -  (l.  -  x3  *  (bb  -  dx))  *  (1.  +  x3  *  bb)**yy 
if  (f3  *  fl  .It.  0.)  then 

X2  =  x3 
else 

xl  =  x3 
endif 

if  (ii.  It.  5)  go  to  9 
xinit  =  0.5  *  (xl  +  x2) 


c 


c       THE  GRID  PARAMETER  IS  DETERMINED  BY  NEWTON-RAPHSON  ITERATION 


c 


xl  =  xinit 
iter  =  0 
10   continue 

iter  =  iter  +  1 
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f  =  (i.  +  xl  *  (bb  -  dx))  *  (1.  -  xl  *  bb)**yy 

&  -  (l.  -  xl  *  (bb  -  dx))  *  (1.  +  xl  *  bb)**yy 

df  =  (bb  -  dx)  *  (1.  -  xl  *  bb)**yy 

&    +  (bb  -  dx)  *  (1.  +  xl  *  bb)**yy 


bb  -  dr)  *  (1.  +  xl  *  (bb  -  dx))  *  (1.  -  xl  *  bb)**yl 
(bb  -  dr)  *  (1.  -  xl  *  (bb  -  dx))  *  (1.  +  xl  *  bb)**yl 


& 
& 

x2  =  xl  -  f  /  df 

xeps  =  abs  (xl  -  x2) 

xl  =  x2 

if  (xeps  .gt.  epps  .and.  iter  .It.  10)  goto  10 

zz  =  xl  *  bb 

hk  =  0.5  /  zz  *  log((l.  +  zz)  /  (1.  -  zz)  ) 

subroutine  chiptemp 
C      THIS  SUBROUTINE  CALCULATES  THE  CHIP  TEMPERATURES  (DIMENSIONAL) 

c 

implicit  real*8  (a-h,o-z) 

common/ng/dxx(40) ,dyy(40) ,dzz(40) ,x(40) ,y(40) ,z(40)  , 
&  dxxs(40) ,dyys(40) ,dzzs(40) 

common/parm/ra , pr , sorsum , dt  ime , xper , rol 1 , dt_inv , xtime 

common / to 1/small , eps, sormax 

common/dims/h , wth , bth , hchip , wchip , bchip , bsub 

common/count/nt,nmax, imax, itmax,  krun,  npnnt 

common/mscn/smp(12, 16, 16) ,resorm(93) , 
&  du(12,16,16) ,dv(12,16,16) , 

&  dw(12, 16,16) ,pp(12, 16, 16) 

common/ limits/ni , nipl , niml , nj , njpl , njml , nk, nkpl , nkml , 
&  nip2,njp2,nkp2, iter,nnmax 

common/fv_init/tod(12,16,16) , uod ( 12 , 16 , 16) , vod ( 12 , 16 , 16) , 
&  wod(12,16, 16) ,pod(12, 16, 16) 

common/ fv_cur/t( 12, 16,16) , u ( 12 , 16 , 16) , v ( 12 , 16 , 16) , 
&  w(12,16,16) ,p(12,16,16) 

COMMON /RHOCP/RHS,RHC 

common/condu/alpha(12,16,16) , rho ( 12 , 16 , 16) , visco ( 12  16 , 16) 
common/chip/ibgn,iend, jbgn(3) , jend(3) ,kbgn(3) ,kend(3) ,nchp, 
&  isub, ichip, jchip,kchip 

common/diff/alc,als,thot,tcool,tavg 
common/array/njchip, nkchip 
COMMON/POWER/QQQ,QCOND 

C 

dimension  tt (3 , 3 ) , vol ( 3 , 3 ) 

c 
C 

DO  15  M=1,NJCHIP 

DO  15  N=l, NKCHIP 
VOL(M,N)  =  0. 
TT(M,N)  =  0. 
15   CONTINUE 


C 


do    20   m=l,njchip 
do    20    n=l, nkchip 

do    10    i=ibgn,iend 

do    10    j=jbgn(m) , jend(m) 
do    10   k=kbgn(n) ,kend(n) 

tt(m,n)    =   tt(m,n)    +   t(i,j,k)    *   dxx(i) 
i  *    dyy(j)    *    dzz(k) 

vol(m,n)    =   vol(m,n)    +   dxx(i)    *    dyy(j) 
r  *    dzz(k) 
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10        continue 

TT(M,N)  =  TT(M,N)  /  VOL(M,N)  *  QCOND 
20   continue 

print  *,'THE  TEMPERATURE  RISES  IN  THE  CHIP  ARE:' 
WRITE(6,200)  TT 
200   FORMAT  (4(F17.6)) 
end 
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